Handling Spatial Projections & CRS
Last updated on 2025-02-18 | Edit this page
Estimated time: 30 minutes
Overview
Questions
- What do I do when vector data do not line up?
Objectives
After completing this episode, participants should be able to…
- Plot vector objects with different CRSs in the same plot.
Working with spatial data from different sources
In this episode, we will work with a different dataset containing Dutch municipal boundaries. We start by reading the data and plotting it.
R
municipal_boundary_NL <- st_read("data/nl-gemeenten.shp")
OUTPUT
Reading layer `nl-gemeenten' from data source
`/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/nl-gemeenten.shp'
using driver `ESRI Shapefile'
Simple feature collection with 344 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
Projected CRS: Amersfoort / RD New
R
ggplot() +
geom_sf(data = municipal_boundary_NL) +
labs(title = "Map of Contiguous NL Municipal Boundaries") +
coord_sf(datum = st_crs(28992))

We can add a country boundary layer to make it look nicer. If we specify a thicker line width using size = 2 for the country boundary layer, it will make our map pop! We read the country boundary from a different file.
R
country_boundary_NL <- st_read("data/nl-boundary.shp")
OUTPUT
Reading layer `nl-boundary' from data source
`/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/nl-boundary.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
Projected CRS: Amersfoort / RD New
R
ggplot() +
geom_sf(data = country_boundary_NL,
color = "gray18",
linewidth = 2) +
geom_sf(data = municipal_boundary_NL,
color = "gray40") +
labs(title = "Map of Contiguous NL Municipal Boundaries") +
coord_sf(datum = st_crs(28992))

We confirm that the CRS of both boundaries is 28992.
R
st_crs(municipal_boundary_NL)$epsg
OUTPUT
[1] 28992
R
st_crs(country_boundary_NL)$epsg
OUTPUT
[1] 28992
We read the municipal boundary of Delft and make sure that it is in the same CRS as the country-level municipal boundaries and country boundary layers.
R
boundary_Delft <- st_read("data/delft-boundary.shp")
OUTPUT
Reading layer `delft-boundary' from data source
`/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/delft-boundary.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 4.320218 ymin: 51.96632 xmax: 4.407911 ymax: 52.0326
Geodetic CRS: WGS 84
R
st_crs(boundary_Delft)$epsg
OUTPUT
[1] 4326
R
boundary_Delft <- st_transform(boundary_Delft, 28992)
R
ggplot() +
geom_sf(data = country_boundary_NL,
linewidth = 2,
color = "gray18") +
geom_sf(data = municipal_boundary_NL,
color = "gray40") +
geom_sf(data = boundary_Delft,
color = "purple",
fill = "purple") +
labs(title = "Map of Contiguous NL Municipal Boundaries") +
coord_sf(datum = st_crs(28992))

Challenge: Plot multiple layers of spatial data
Create a map of South Holland as follows:
- Import
nl-gemeenten.shp
and filter only the municipalities in South Holland. - Plot it and adjust line width as necessary.
- Layer the boundary of Delft onto the plot.
- Add a title.
- Add a legend that shows both the municipal boundaries (as a line) and the boundary of Delft (as a filled polygon).
R
boundary_ZH <- municipal_boundary_NL %>%
filter(ligtInPr_1 == "Zuid-Holland")
R
ggplot() +
geom_sf(data = boundary_ZH,
aes(color ="color"),
show.legend = "line") +
scale_color_manual(name = "",
labels = "Municipal Boundaries in South Holland",
values = c("color" = "gray18")) +
geom_sf(data = boundary_Delft,
aes(shape = "shape"),
color = "purple",
fill = "purple") +
scale_shape_manual(name = "",
labels = "Municipality of Delft",
values = c("shape" = 19)
) +
labs(title = "Delft location") +
theme(legend.background = element_rect(color = NA)) +
coord_sf(datum = st_crs(28992))

Projecting layers
Note that ggplot2
may reproject the layers on the fly
for visualisation purposes, but for geoprocessing purposes, you still
need to reproject the layers explicitly with
st_transform()
. This will become clear in a later episode when we perform GIS
operations.
Export a shapefile
To save a file, use the st_write()
function from the
sf
package. Although sf
guesses the driver
needed for a specified output file name from its extension, this can be
made explicitly via the driver
argument. In our case
driver = "ESRI Shapefile"
ensures that the output is
correctly saved as a .shp
file.
R
st_write(leisure_locations_selection,
"data/leisure_locations_selection.shp",
driver = "ESRI Shapefile"
)
Key Points
-
ggplot2
automatically converts all objects in a plot to the same CRS. - For geoprocessing purposes, you still need to reproject the layers you use to the same CRS.
- You can export an
sf
object to a shapefile withst_write()
.