Spatial R: Using the sf package

Working with spatial data

Working with spatial data in R is easy! In this article, we’ll be using the Simple Features for R package, also known as the sf package which makes it straightforward to manipulate and analyze spatial data. If you do not have the sf or any another package, you can install it using the install.packages() function. After installing these packages, you can load them for each R session using library() as we will do here. Let’s go ahead and load the various R packages we’ll use in this article.

library(sf)          #For working with spatial data
library(tidyverse)   #For data wrangling and ggplot2
library(ggspatial)   #For spatial data with ggplot2
library(tigris)      #For downloading spatial data directly into R

You love cities and recently moved to Charlottesville, Virginia, for school but you’re excited to explore all the places Virginia has to offer. However, you can’t find a map that provides all the information you want, so you need to use spatial information to build your own map to help you navigate around. Let’s start by reading in our first layer, the boundary or shape of Virginia. This layer is a polygon and it is contained in a shapefile. A shapefile is simply a format type for storing geographic data and has a .shp extension. A shapefile is part of a network of multiple files that all work together, so it’s important to store them in the same folder. Click here to download the Virginia spatial files used in this article if you would like to follow along.

We will read in a shapefile of a Virginia polygon using st_read() and check the geometry type by calling st_geometry_type() on our spatial data. Checking the geometry will tell us the type of spatial data we are working with (ex. point, linestring). The polygon we will use for Virginia is from the Virginia Department of Emergency Management.

#Read in Virginia polygon 
virginia <- st_read("VirginiaBoundary.shp") 

#Check geometry type 
as.character(st_geometry_type(virginia))[1]

[1] "POLYGON"

We can take a look at our Virginia boundary by calling the plot() function on our spatial data. We need to tell R to specifically plot the geometry column that contains the spatial data, otherwise, R will plot the Virginia boundary for every column in our data frame.

plot(virginia$geometry)

Great! We have the outline of Virginia. Let’s add our city locations by making a data frame with the cities that we are interested in visiting. We will take a look at the structure of this data frame using str().

#Data frame of cities 
cities <- data.frame(
  City = c("Richmond", "Virginia Beach", "Norfolk", "Chesapeake", "Arlington", "Charlottesville"),
  Latitude = c(37.5413, 36.8529, 36.8508, 36.7682, 38.8797, 38.0293),
  Longitude = c(-77.4348, -75.9780, -76.2859, -76.2875, -77.1080, -78.4767)) 

#Examine structure of data frame
str(cities)

'data.frame':   6 obs. of  3 variables:
 $ City     : chr  "Richmond" "Virginia Beach" "Norfolk" "Chesapeake" ...
 $ Latitude : num  37.5 36.9 36.9 36.8 38.9 ...
 $ Longitude: num  -77.4 -76 -76.3 -76.3 -77.1 ...

We can see that our data frame has 6 observations of 3 variables or columns. In our data frame, we have a column for the city name, as well as the latitude and longitude of the city. The latitude and longitude columns define our coordinates and provide locations of these places on earth. Notably, our cities are characters while our latitude and longitude are numeric values. But, how does R know these numbers should represent locations? To tell R this data is spatial information, we are going to convert these values to a column that contains a geometry. Like our polygon of Virginia, our geometry column contains the data structure to represent spatial features. To do this conversion, we will use the function st_as_sf() from the sf package. We will input the data frame to use, the columns to convert to a geometry (our Latitude and Longitude columns) and also tell R to not remove the original Longitude and Latitude columns by setting the argument remove = FALSE. You can choose to remove the numeric Longitude and Latitude columns from the data frame by leaving this argument out or setting it to TRUE, which is the default. We will also start by setting crs = 4326, a common global system to reference our locations. We will review the crs argument in the next section.

#Convert to sf object
cities.sf <- 
  st_as_sf(cities, coords = c("Longitude", "Latitude"), remove = FALSE, crs = 4326)

#Plot the cities
plot(cities.sf$geometry)

When we work with data frames containing latitude and longitude information, we have the opportunity to understand the spatial context of the data. Keep in mind that these latitude and longitude coordinates are point data representing specific locations on earth. However, looking at the city points alone doesn’t give us the spatial context we need, so lets add the Virginia boundary as a layer. Before we can put these 2 layers together, we need to look at how our data is represented because we can measure spatial data in different ways.

We will check this by calling st_crs() from the sf package to get the structure of our spatial data. This function will give us information on the coordinate reference system (CRS), which we will explore in the next section.

st_crs(virginia)$epsg

[1] 4269

st_crs(cities.sf)$epsg

[1] 4326

Essentially, the way we have represented our layers are different. What does this mean and how can we fix it?

Defining spatial data

Imagine you are now working in a coffee shop and the manager asks you to measure the bar length so they can order a new one since this bar has gotten old. You tell your manager it is four. 4 feet? 4 meters? Certainly it can’t be very small coming in at 4 cm or very large at 4 miles! If it is truly 4 feet, we can say it is also 1.2 meters (4 ft = 1.2 m) while if it is actually 4 meters, the bar is longer at 13 feet and can probably fit a few more people. Both of these options seem reasonable to the manager since they have not been into the shop before, but which is it? It certainly can’t be both since the length of the bar hasn’t changed, but rather the unit of measuring the bar has.

Similarly, we have different ways to present the spatial information associated with our data. We will explore ways to show our data by going over two important components of spatial data: (a) datums and (b) projections. We’ll also cover how to represent them using (a) coordinate reference systems and (b) EPSG codes.

Remember that coordinates are a pair of numbers that tell us the location of something. A datum specifies the frame for our coordinates. This means that a datum provides the foundation for defining where things are relative to one another. In other words, it tells us the “origin and orientation of latitude and longitude lines” (esri, Datums). Coordinates are associated with units (ex. meters, degrees) and are also different for different datums, so defining the correct datum is important!

Phew! We can represent locations on earth now and have identified a means to do so. But, we still face a challenge - the earth isn’t flat. It’s a sphere which is 3-dimensional while we want to create a map that is only 2-dimensional. We can manage this problem by setting a projection. Defining a projection allows us to represent something that is 3-dimensional on a 2-dimensional plane. Essentially, a projection tells us how the 3-dimensional space will be projected, or drawn, onto a 2-dimensional surface.

So, the datum gives us a reference frame for the location of things while the projection tells us how to draw them on a flat surface. We can tie these components together by using a CRS (coordinate reference system). A CRS offers a standardized method to define locations and can be referenced by an EPSG code which is a unique identifier for geographic objects. You can find more on EPSG codes by visiting this website. We need our CRS to be consistent across all layers so that they align properly. If the CRS is different, you will likely encounter an error in R, indicating that the CRS of layer1 doesn’t match the CRS of layer2.

Now that we know about the structure of our spatial data, we need to ensure our Virginia polygon and city points have the same CRS. We will do this by calling st_transform() on our data. This sf function will take our spatial data and transform it to the CRS we set. We will use the EPSG 26918 which is a projected coordinate system commonly used in North America. We will then check that our layers have the same CRS by seeing if the CRS of virginia.crs is equal to (==) the CRS of cities.crs.

#Transform CRS
virginia.crs <- st_transform(virginia, crs = 26918)
cities.crs <- st_transform(cities.sf, crs = 26918)

#Check that the layers have the same CRS
st_crs(virginia.crs) == st_crs(cities.crs)

[1] TRUE

If you run into further issues or want to explore the data more, you can also troubleshoot by checking the coordinate values by calling head(st_coordinates(df)) on your data frame.

Plotting with ggplot

Now that we’ve brought in our data, converted it to a spatial object, and made sure our layers align, let’s create our map. We’ll use ggplot() and call geom_sf() which tells R that we’re working with spatial data. When making a map, it’s also crucial to consider the order of our layers. If we place our cities data (cities.crs) before the Virginia polygon (virginia.crs), the cities will end up under the Virginia polygon and we won’t be able to see them.

ggplot() +
  #Add Virginia polygon
  geom_sf(data = virginia.crs) +
  #Add cities layer
  geom_sf(data = cities.crs, aes(color = City, label = City), size = 3) +
  #Add titles
  labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
  #Add theme
  theme_bw()

Since the cities we want to visit are in the eastern region of the state, we’ll ‘zoom in’ on this area by defining a bounding box around our cities. A bounding box defines a rectangle that encompasses the extent of our features of interest. We will get our bounding box of the cities extent and crop virginia.crs to this region by calling st_crop().

#Get bounding box of the cities 
cities.extent <- st_bbox(cities.crs)

Additionally, we’ll add some custom elements to our map, like a north arrow and a scale bar, using the ggspatial R package.

#Plot using the packages ggplot and ggspatial 
ggplot() +
  #Crop Virginia boundary to spatial extent of cities and add Virginia layer
  geom_sf(data = st_crop(virginia.crs, cities.extent)) +
  #Add cities layer
  geom_sf(data = cities.crs, aes(color = City, label = City), size = 3) +
  #Add scale bar to bottom left from ggspatial
  annotation_scale(location = "bl", height = unit(.25, "cm"), 
                   width = unit(1, "cm"), pad_x = unit(0.3, "in"), 
                   pad_y = unit(0.3, "in")) +
  #Add north arrow to bottom left from ggspatial
  annotation_north_arrow(height = unit(1, "cm"), width = unit(1, "cm"),
                         which_north = "true", location = "bl", 
                         pad_x = unit(0.4, "in"), pad_y = unit(0.5, "in")) +
  #Add titles
  labs(x = "Longitude", y = "Latitude", title = "Virginia Cities") +
  #Add theme
  theme_bw()

Extracting information in spatial data

There is still a lot of information we can gather from this data outside of the map itself! How far from the Virginia border are the cities? What is the distance from each city to Charlottesville? What is the area of each city that we can explore? Answering these questions is easy using sf! We will first get the distance from each city to the border using the st_distance() function on the boundary of the Virginia polygon (st_boundary()). The st_distance() function calculates the distance between two spatial objects. Here, we calculate the Euclidean distance, or the shortest distance between two spatial objects.

#Distance from each city to the border of Virginia 
dist.border <-  st_distance(cities.crs, st_boundary(virginia.crs))

#See results
print(dist.border)

Units: [m]
           [,1]
[1,] 85469.5113
[2,]   121.7268
[3,] 11906.0329
[4,] 20091.7204
[5,]  3054.4664
[6,] 82565.9571

When we take a look at the distances, we see [m] at the top, indicating that these results are in the units of meters. If you would like to convert these values to a numeric class, you can do so using as.numeric() on the distance column, or you can add them to your data frame for later use by calling cities.crs$distance_to_border <- dist.border. From our results, we can see that Virginia Beach is closest to the Virginia border while Richmond is farthest.

However, we are in Charlottesville, and so we want how far each city is from Charlottesville. We calculate these distances and add them to our cities.crs data frame.

#Calculate distance from Charlottesville to each city 
dist.cville <- 
  st_distance(cities.crs, filter(cities.crs, City == "Charlottesville")) %>% 
  cbind(cities.crs)

#Look at first 3 distances
head(dist.cville, n = 3)

Simple feature collection with 3 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 284889.6 ymin: 4079001 xmax: 412812.9 ymax: 4157709
Projected CRS: NAD83 / UTM zone 18N
             .           City Latitude Longitude                 geometry
1 106615.8 [m]       Richmond  37.5413  -77.4348 POINT (284889.6 4157709)
2 256808.3 [m] Virginia Beach  36.8529  -75.9780 POINT (412812.9 4079001)
3 233902.9 [m]        Norfolk  36.8508  -76.2859 POINT (385359.8 4079093)

Since we’ve figured out how far these cities are from us, let’s now look at some navigation routes and the area of each city to estimate the size and the space available for us to explore. We will start by obtaining county data for Virginia using the tigris package, which provides access to geographic and spatial information from the U.S. Census Bureau. We’ll retrieve layers for the county boundaries (counties()) and primary roads (primary_roads()). Our road layer is a linestring or a sequence of continuous lines.

Keep in mind that it’s important to check the CRS of each layer and transform them if they do not match; otherwise, we’ll encounter errors when plotting the layers due to mismatched projections.

#Get county data for Virginia from tigris package and transform the CRS 
va.counties <- counties(state = "VA") %>% 
  st_transform(crs = 26918)

#Get road data and transform the CRS
roads <- primary_secondary_roads(state = "VA", year = 2021) %>% st_transform(crs = 26918)

#Plot with ggplot 
ggplot() +
  #Crop Virginia counties to cities extent and plot 
  geom_sf(data = st_crop(va.counties, cities.extent), color = "gray100", fill = "gray88", size = .25) +
  #Crop roads to cities extent and plot 
  geom_sf(data = st_crop(roads, cities.extent), color = "black", size = .3) +
  #Plot cities 
  geom_sf(data = cities.crs, aes(color = City), size = 4) +
  #Add titles 
  labs(x = "Longitude", y = "Latitude", title = "Map of Virginia Cities") +
  #Adjust theme
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

Now that we can visualize the boundary of each city and potential navigation routes to get there, we can calculate the area of our cities using st_area(). We will calculate the area for each of our 6 cities and create a column, area using mutate() to hold the calculated areas.

#Calculate area of each city we want to visit 
county.area <- va.counties %>% 
  select(NAME, geometry) %>%         #Select NAME and geometry columns
  filter(NAME %in% cities$City) %>%  #Filter to rows where NAME is in the City column of the cities data frame
  mutate(area = st_area(geometry))   #Create a column of areas using st_area()

#Output county name and area 
county.area %>%
  select(NAME, area)

Simple feature collection with 7 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 190674 ymin: 4045325 xmax: 428626.1 ymax: 4311631
Projected CRS: NAD83 / UTM zone 18N
             NAME             area                       geometry
1 Charlottesville   26620749 [m^2] MULTIPOLYGON (((190676.7 42...
2        Richmond  162115691 [m^2] MULTIPOLYGON (((270200.2 41...
3      Chesapeake  908525938 [m^2] MULTIPOLYGON (((366512.5 40...
4       Arlington   67607783 [m^2] MULTIPOLYGON (((311606.8 43...
5         Norfolk  249559097 [m^2] MULTIPOLYGON (((380155.2 40...
6  Virginia Beach 1287748687 [m^2] MULTIPOLYGON (((390519.1 40...
7        Richmond  560255711 [m^2] MULTIPOLYGON (((358946.8 42...

Here, we can see the area of each city in \(m^2\). We can also see that the geometry type of our counties is a multipolygon, or a collection of polygons. This makes sense since each county that we have is a polygon shape, and we have multiple of them! We can also see the coordinates of our bounding box in Bounding box as well the CRS used. Remember, we can also get the spatial information by calling st_crs(df)$datum for the datum or st_crs(df)$proj for the projection.

st_crs(county.area)$datum

[1] "NAD83"

st_crs(county.area)$proj

[1] "utm"

We could easily convert our area values to square miles or square kilometers by applying the correct conversion factor, if needed.

Conclusion

We’ve explored different types of spatial data, ensured proper formatting, and covered several techniques for visualization and data analysis. However, the ability to use spatial information does not end here! You can use distance or area as a predictor variable in a model, visualize income by city on a map, or simply plot your latitude and longitude points for a simple representation of your data locations. You can join datasets based on geographic features (st_join()), find where features intersect (st_intersects()), or even do geospatial statistics. There are also opportunities to use many layers or data sources and include different spatial data types, like we’ve done here by including point, linestring, and polygon data.

To further explore the sf package, check out this R spatial github, and if you’re interested in creating interactive maps, check out our article Data Scientist as Cartographer: An Introduction to Making Interactive Maps in R with Leaflet!

References

• Dunnington, D. (2023). ggspatial: Spatial Data Framework for ggplot2. R package version 1.1.9, https://CRAN.R-project.org/package=ggspatial.
• Esri https://desktop.arcgis.com/en/arcmap/latest/map/projections/datums.htm
• Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016
• R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Virginia Department of Emergency Management (2024). https://vgin.vdem.virginia.gov/datasets/777890ecdb634d18a02eec604db522c6/about.
• Walker, K. (2023). tigris: Load Census TIGER/Line Shapefiles. R package version 2.0.4, https://CRAN.R-project.org/package=tigris.
• Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, LD., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, TL., Miller, E., Bache, SM., Müller, K., Ooms, J., Robinson, D., Seidel, DP., Spinu, V., Takahashi, K., …Yutani, H. (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.