Handling Spatial Projections & CRS
Last updated on 2024-12-03 | Edit this page
Overview
Questions
- What do I do when vector data don’t 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
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!
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))
R
# st_crs(point_Delft)
R
st_crs(municipal_boundary_NL)$epsg
OUTPUT
[1] 28992
R
st_crs(country_boundary_NL)$epsg
OUTPUT
[1] 28992
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 province 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))
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. - Still be aware of the CRS and extent for each object.
- You can export an
sf
object to a shapefile withst_write()
.