Key Learnings from, and Solutions to the exercises in Chapter 7 of the book Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow.
Geocomputation with R
Textbook Solutions
Author
Aditya Dahiya
Published
January 10, 2025
Code
library(sf) # Simple Features in Rlibrary(terra) # Handling rasters in Rlibrary(tidyterra) # For plotting rasters in ggplot2library(magrittr) # Using pipes with raster objectslibrary(tidyverse) # All things tidy; Data Wranglinglibrary(spData) # Spatial Datasetslibrary(patchwork) # Composing plotsgreenland <- rnaturalearth::ne_countries(geounit ="Greenland",scale ="medium",returnclass ="sf" ) |>select(admin, name, iso_a3, geometry)sysfonts::font_add_google("Saira Condensed", "body_font")showtext::showtext_auto()theme_set(theme_minimal(base_family ="body_font",base_size =14,base_line_size =0.2 ) +theme(text =element_text(colour ="grey20" ),plot.title =element_text(size =24,margin =margin(0.1, 0.1, 0.1, 0.1, "lines") ),plot.subtitle =element_text(size =16,margin =margin(0.1, 0.1, 0.1, 0.1, "lines"),lineheight =0.3 ) ))
Note
For the first few sections, I will use the map of Greenland as an example, since it is an extreme case (far away from equator, so a good use for CRS projections examples), as shown in Figure 1.
7.1 Introduction
Coordinate Reference Systems (CRSs): Two main types:
Geographic CRS: A reference system that defines locations on the Earth using a spherical or ellipsoidal model. It specifies positions in terms of latitude and longitude, with angular measurements relative to the Earth’s center. Uses longitude/latitude (units in degrees)
Projected CRS: A reference system that transforms the Earth’s curved surface into a flat map. It uses linear units like meters and focuses on minimizing distortions in specific aspects such as distance, area, or shape, depending on the map’s purpose. Uses meters from a datum (projected system).
Aspect
Geographic CRS
Projected CRS
Definition
Represents locations on the Earth’s surface using a spherical or ellipsoidal model.
Represents locations on a flat, two-dimensional surface using a projected system.
Units of Measurement
Uses angular units (degrees of longitude and latitude).
Uses linear units (meters, feet, etc.).
Coordinate Axes
Longitude (x-axis) and Latitude (y-axis).
X (east-west) and Y (north-south) in a flat plane.
Purpose
Ideal for global-scale mapping and geodetic calculations.
Ideal for local or regional mapping, where accurate distances, angles, and areas are required.
Accuracy
Retains accurate angular relationships but distorts distances and areas.
Distorts angular relationships but preserves distances, areas, or shapes (depending on projection).
Examples
WGS84 (used in GPS), NAD83.
UTM (Universal Transverse Mercator), State Plane.
Key Feature
Geocentric, based on Earth’s shape as an ellipsoid.
Flat, derived by projecting the Earth onto a 2D surface using a mathematical formula.
Applications
Satellite data, GPS, global navigation systems.
Engineering projects, land use planning, local/regional maps.
Focus of the Chapter:
Setting and transforming CRSs in geographic data.
Highlighting issues caused by ignoring CRSs, especially with lon/lat data.
Practical Relevance:
Many projects don’t require conversion between CRSs but knowing whether your data is in a projected or geographic CRS is critical.
Understanding CRSs prevents errors and ensures smooth handling of geographic data.
CRS in R spatial objects: Coordinate Reference Systems (CRS) are queried and set in vector and raster geographic data using functions from the sf(Pebesma and Bivand 2023a) and terra(Hijmans 2024a) packages.
Vector data CRS handling:
Use st_crs() to query CRS in vector objects (e.g., EPSG:4326 for WGS 84).
CRS metadata includes User input (e.g., EPSG code, proj-string) and wkt (WKT string with complete CRS details).
Additional CRS details can be retrieved: The st_crs() function from the sf package in R provides detailed information about the Coordinate Reference System (CRS) of spatial objects.
Attribute
Description
Example Usage
IsGeographic
Indicates if the CRS is geographic (TRUE) or projected (FALSE).
st_crs(new_vector)$IsGeographic
units_gdal
Specifies the units of measurement used in the CRS, as recognized by GDAL.
st_crs(new_vector)$units_gdal
srid
Provides the Spatial Reference System Identifier (SRID) associated with the CRS, if available.
st_crs(new_vector)$srid
proj4string
Returns the PROJ.4 string representation of the CRS, offering a concise textual description.
st_crs(new_vector)$proj4string
WktPretty
Provides a formatted Well-Known Text (WKT) representation of the CRS for easier readability.
st_crs(new_vector)$WktPretty
InvFlattening
Provides the inverse flattening value of the ellipsoid, indicating its compression.
st_crs(new_vector)$InvFlattening
Name
Retrieves the official name of the CRS.
st_crs(new_vector)$Name
epsg
Returns the EPSG code associated with the CRS, if available.
st_crs(new_vector)$epsg
yx
Indicates whether the axis order is latitude-longitude (TRUE) or longitude-latitude (FALSE).
st_crs(new_vector)$yx
ud_unit
Returns the units object associated with the CRS, or NULL if units are missing.
st_crs(new_vector)$ud_unit
axes
Provides information about the axes of the CRS, including their names and units.
st_crs(new_vector)$axes
Code
# The WGS 84 / EPSG 4326 projection (default projection in {rnaturalearth})st_crs("EPSG:4326")$units_gdal## [1] "degree"st_crs("EPSG:4326")$IsGeographic## [1] TRUEst_crs("EPSG:4326")$proj4string## [1] "+proj=longlat +datum=WGS84 +no_defs"st_crs("EPSG:4326")$axes## name orientation## 1 Geodetic latitude 1## 2 Geodetic longitude 3# The EPSG:3857 Web Mercator Projection (used for maps and distances)st_crs("EPSG:3857")$units_gdal## [1] "metre"st_crs("EPSG:3857")$IsGeographic## [1] FALSEst_crs("EPSG:3857")$proj4string## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"st_crs("EPSG:3857")$axes## name orientation## 1 Easting 3## 2 Northing 1# Robinson Projection ESRI:54030 st_crs("ESRI:54030")$units_gdal## [1] "metre"st_crs("ESRI:54030")$IsGeographic## [1] FALSEst_crs("ESRI:54030")$proj4string## [1] "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"st_crs("ESRI:54030")$axes## name orientation## 1 Easting 3## 2 Northing 1
Use st_set_crs() to manually set or correct a CRS. Both st_set_crs() and the assignment method st_crs(object) <- value are used to define or modify the Coordinate Reference System (CRS) of spatial objects.
st_set_crs() is a function call, which can be advantageous when using pipes (|>) in a workflow.
Both methods assign CRS metadata to the spatial object without altering the actual coordinate values or geometries.
To transform the coordinates to a different CRS, the st_transform() function should be used.
Assigning a CRS is akin to labeling data with its existing coordinate system, whereas transforming reprojects the data into a new coordinate system.
Explaining the difference between the st_set_crs() and st_transform(). The figures show (a) The default map of Greenland as obtained from {rnaturalearth} with the CRS WGS84 (i.e. EPSG:4326). (b) The transformation using st_transform() which produces correct transformation. (c) The incorrect method using st_set_crs() which simply retains all data, but assumes them to be in the new assigned CRS. This leads to wrong latitude and longitude numbers. Thus, the only use of st_set_crs() is for assigning a CRS to a raw data, which we know for sure from outside sources.
Validate CRS with st_is_longlat().
The function st_is_longlat() in the sf package is used to check whether the Coordinate Reference System (CRS) of a spatial object is based on a geographic coordinate system (latitude and longitude). It’s Return Value can be: —
TRUE: The object has a geographic CRS (e.g., EPSG:4326).
FALSE: The object has a projected CRS.
NA: The CRS is undefined (e.g., no metadata about the coordinate system is set).
Raster data CRS handling:
Use crs() from the {terra} package to query and set CRS for raster objects.
Outputs CRS as a WKT string.
Set CRS with EPSG codes, WKT strings, or CRSs from other objects (e.g., crs(my_raster) = "EPSG:26912").
Important considerations:
st_crs() and crs() only set metadata; they do not modify coordinate values or geometries.
Missing CRS in datasets (e.g., NA output in st_is_longlat()) indicates that CRS must be explicitly defined (e.g., using st_set_crs()).
CRS assumptions: R avoids assumptions about CRS (e.g., unlike GeoJSON which defaults to EPSG:4326).
7.4 Geometry Operations on Projected and Unprojected Data
GEOS is used for projected CRS or CRS-absent data.
Default Behavior: S2 is enabled by default for geographic data but can be disabled using sf::sf_use_s2(FALSE).
CRS Importance in Buffers: Buffers created on unprojected data (e.g., lon/lat) can result in distorted outputs (for example A.1 in Figure 2) unless spherical geometry (S2) is used. Better still, always use a projected CRS for accurate distance / buffers calculation (as C.1 in Figure 2).
Key Learnings:
Incorrect buffer creation: Buffers on geographic data (Long / Lat data) without a specified CRS (thus, using GEOS) result in outputs in degrees (as A.1 & A.3 in Figure 2). Note: Lon/lat buffers without S2 or projection may appear elongated (e.g., along the north-south axis).
Correct buffer creation: S2 buffers (with Geographic CRS, i.e. Long / Lat data) calculate accurate distances (as B.1 in Figure 2) but may create jagged boundaries (as B.3 in Figure 2) (this can be tuned with max_cells or nQuadSegs).
Ideal buffer creation: Use projected CRSs (as in C.1, C.2 & C.3 in Figure 2)) like British National Grid (EPSG:27700) for consistent distance-based operations, for example, in metres (m).
Note
Note: S2 buffers improve results but remain less precise than using projected CRSs and GEOS.
Visual Outcomes:
Buffers based on S2 and projected CRSs produce more accurate, equidistant boundaries compared to lon/lat inputs without S2.
S2-derived buffers differ from projected GEOS buffers, and resolution depends on max_cells (default: 1000). (as B.3 in Figure 2)
Best Practices:
Use sf::st_crs() to check the CRS and verify units (e.g., degrees vs. meters) of the CRS. Reproject data to a projected CRS using suitable EPSG codes for precise geometric operations.
For distance-based operations like buffering, reproject data to a projected CRS.
Adjust max_cells in S2 for performance versus resolution trade-offs. (as B.1 in Figure 2)
Figure 2 shows an example (buffer zone around London)
Code
overall_title <-"When and How to Use S2 and GEOS Engines for Spatial Operations"overall_subtitle <-"GEOS handles projected or CRS-absent data; S2 defaults for geographic data but is optional."bts <-16# Base Text Sizetheme_set(theme_minimal(base_family ="body_font",base_size = bts,base_line_size =0.2 ) +theme(text =element_text(colour ="grey20" ),plot.title =element_text(size =1.9* bts,margin =margin(0.2, 0.1, 0.1, 0.1, "lines") ),plot.subtitle =element_text(size =1.3* bts,margin =margin(0.1, 0.1, 0.1, 0.1, "lines") ),plot.margin =margin(0,0,0,0, "mm"),axis.text =element_text(size = bts *0.5 ),axis.ticks.length =unit(0, "mm"),panel.background =element_rect(fill ="#caf0f8",colour ="#caf0f8" ),panel.grid =element_line(linewidth =0.2,colour ="white",linetype =3 ) ))# Create a London data point - the coordiantes in Latitude# and Longitude of London, UKlondon_sans_crs <-data.frame(lon =-0.1, lat =51.5) |>st_as_sf(coords =c("lon", "lat"))# Adding the CRS to London data point: values stay samelondon_with_crs <- london_sans_crs |>st_set_crs("EPSG:4326")# British National Grid# st_crs("EPSG:27700")# Using the British National Grid CRS: values get changedlondon_british_crs <- london_sans_crs |>st_set_crs("EPSG:4326") |>st_transform("EPSG:27700")# Create an approximate bounding box around 200 km around Londonbbox_near_london <- london_sans_crs |>st_set_crs("EPSG:4326") |>st_buffer(dist =200000) |>st_bbox()# Crop and keep coastline map around Londonoutline_map <- rnaturalearth::ne_countries(continent ="Europe",returnclass ="sf",scale ="large" ) |>st_union() |>st_as_sf() |>st_crop(bbox_near_london)base_plot <-ggplot() +geom_sf(data = outline_map,linewidth =0.1,fill =alpha("white", 1) ) +geom_sf(data = london_with_crs,size =3,colour ="#fb8500" ) +coord_sf(expand =FALSE)# Without any CRS -------------------------------------------------row_1_title <-"London Data, without any CRS"row_1_subtitle <-"Data in Long / Lat, but without a specified CRS, leads {sf} to use GEOS engine, as in (a) and (c)"g1 <- base_plot +geom_sf(data =st_buffer(london_sans_crs, dist =1) |>st_set_crs("EPSG:4326"),fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Incorrect: Using GEOS on Long/Lat data",subtitle ="st_buffer(london_sans_crs, dist = 1)" )g2 <- base_plot +geom_sf(data = london_sans_crs |>st_set_crs("EPSG:4326") |>st_buffer(dist =1),fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Incorrect: Using S2 but wrong Buffer units",subtitle ="london_sans_crs |> st_set_crs(`EPSG:4326`) |> st_buffer(dist = 1)" )g3 <- base_plot +geom_sf(data =st_buffer(london_sans_crs, dist =1,nQuadSegs =2) |>st_set_crs("EPSG:4326"),fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Incorrect: Using GEOS with nQuadSegs",subtitle ="st_buffer(london_sans_crs, dist = 1, nQuadSegs = 3)" )row_2_title <-"London with a Geographic CRS `EPSG:4326`"row_2_subtitle <-"Data in Long/Lat, but with CRS speficied, leads {sf}to use S2 (spherical geometry). Distance is in metres, now."sf_use_s2(TRUE)g4 <- base_plot +geom_sf(data =st_buffer(london_with_crs, dist =1e5,max_cells =300),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Correct: Buffer of 100 km using S2",subtitle ="Using S2 engine (spherical geometry); Notice jagged buffer line." )sf_use_s2(FALSE)g5 <- base_plot +geom_sf(data =st_buffer(london_with_crs, dist =1,max_cells =500),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Incorrect: Buffer of 1 degree with GEOS",subtitle ="sf_use_s2(FALSE): Using GEOS engine on geographic CRS." )sf_use_s2(TRUE)g6 <- base_plot +geom_sf(data =st_buffer(london_with_crs, dist =1e5,max_cells =100),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Correct: Buffer of 100 km using S2",subtitle ="Reducing max_cells = 100 to save computation time." )row_3_title <-"With a Projected CRS: British Grid (EPSG:27700)"row_3_subtitle <-"Using London data in Projected CRS results in ideal outputs. Accurate results with GEOS engine."g7 <- base_plot +geom_sf(data =st_buffer(london_british_crs, dist =1e5),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Ideal: 100 km Buffer on projected CRS",subtitle ="london_british_crs |> st_buffer(dist = 1e5)" )g8 <- base_plot +geom_sf(data =st_buffer(london_british_crs, dist =1e5,nQuadSegs =1),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Distance in 4 directions only",subtitle ="nQuadSegs = 1 (Number of segments per quadrant)." )g9 <- base_plot +geom_sf(data =st_buffer(london_british_crs, dist =1e5,nQuadSegs =2),linewidth =1,fill =alpha("#ffb703", 0.5),colour ="#ffb703" ) +labs(title ="Save Computation Time: 100 km Buffer",subtitle ="nQuadSegs = 2 (Number of segments per quadrant)." )# A custom theme for each Row headingstheme_row <-function(...){theme(plot.title =element_text(margin =margin(5,0,1,0, "mm"),size = bts *3.5,hjust =0.5,face ="bold" ),plot.subtitle =element_text(margin =margin(0,0,1,0, "mm"),size = bts *2.5,hjust =0.5,face ="bold" ), ... )}theme_tag <-function(...){theme(plot.tag =element_text(margin =margin(0,0,0,0, "mm"),hjust =0.5,face ="bold",size = bts *1.5 ),plot.tag.position ="top" )}# Compiling Plots for London data without CRSga <- g1 + g2 + g3 +plot_annotation(title = row_1_title,subtitle = row_1_subtitle,theme =theme_row(),tag_levels ="1",tag_prefix ="A." ) &theme_tag()# Save the 3 plots for 1st Rowggsave(plot = ga,filename = here::here("book_solutions", "images", "temp_chapter7-4-1.png"),height =1050,width =2400,unit ="px",bg ="white")# Compiling Plots for London data with a Geographic CRSgb <- g4 + g5 + g6 +plot_annotation(title = row_2_title,subtitle = row_2_subtitle,theme =theme_row(),tag_levels ="1",tag_prefix ="B." ) &theme_tag()# Save the 3 plots for 2nd Rowggsave(plot = gb,filename = here::here("book_solutions", "images", "temp_chapter7-4-2.png"),height =1050,width =2400,unit ="px",bg ="white")# Compiling Plots for Projected (British Grid) CRSgc <- g7 + g8 + g9 +plot_layout(tag_level ="new" ) +plot_annotation(title = row_3_title,subtitle = row_3_subtitle,theme =theme_row(),tag_levels ="1",tag_prefix ="C." ) &theme_tag()# Save the 3 plots for 3rd Rowggsave(plot = gc,filename = here::here("book_solutions", "images", "temp_chapter7-4-3.png"),height =1050,width =2400,unit ="px",bg ="white")# Compile the imageslibrary(magick)ga <-image_read(here::here("book_solutions", "images", "temp_chapter7-4-1.png"))gb <-image_read(here::here("book_solutions", "images", "temp_chapter7-4-2.png"))gc <-image_read(here::here("book_solutions", "images", "temp_chapter7-4-3.png"))image_append(c(ga, gb, gc), stack =TRUE) |>image_write(path = here::here("book_solutions", "images", "chapter7-4.png"))# Clean up the temporary filesunlink(here::here("book_solutions","images","temp_chapter7-4-1.png"))unlink(here::here("book_solutions","images","temp_chapter7-4-2.png"))unlink(here::here("book_solutions","images","temp_chapter7-4-3.png"))
Figure 2: The different kinds of buffering - incorrect, correct and ideal using {sf} in R relying on two engines - Google’s S2 spherical geometry, or GEOS engine.
Note
This code below calculates the distances between lines of longitude (meridians) and latitude at different geographic points to demonstrate the distortion inherent in geographic coordinate systems (lon/lat CRSs).
# Load the required packagelibrary(geosphere)# 1. Distance between two meridians (longitude lines) at the equator (0° latitude)# From (0°E, 0°N) to (1°E, 0°N)distance_meridian_equator <-distGeo(c(0, 0), c(1, 0)) cat("Distance between meridians at the equator:", round(distance_meridian_equator/1000, 2), "km\n")## Distance between meridians at the equator: 111.32 km# 2. Distance between two meridians at the latitude of London (51.5°N)# From (0°E, 51.5°N) to (1°E, 51.5°N)distance_meridian_london <-distGeo(c(0, 51.5), c(1, 51.5)) cat("Distance between meridians at London's latitude:", round(distance_meridian_london/1e3, 2), "km\n")## Distance between meridians at London's latitude: 69.44 km# 3. Distance between two parallels (latitude lines) at any latitude# From (0°E, 0°N) to (0°E, 1°N)distance_latitude <-distGeo(c(0, 0), c(0, 1)) cat("Distance between latitudes:", round(distance_latitude/1e3, 2), "meters\n")## Distance between latitudes: 110.57 meters
Thus, the distance between meridians varies with latitude, shrinking from approximately 111 km at the equator to zero at the poles, while lines of latitude remain equidistant.
7.5 When to Reproject?
CRS Assignment: In practical daily work, CRSs are typically set automatically when data is read; hoever, if required, manual setting can be done using st_set_crs().
The main task, in daily work, involves transforming objects between CRSs to enable specific tasks. Transformation Necessity:
For projects with web mapping, we need to use EPSG:4326 (a geographic CRS).
For projects needing Spherical Geometry operations, use Geographic CRS.
Tranformation into another Projected CRS
Projects needing planar geometry operations, like creating buffer zones with smooth edges.
Required when comparing or combining objects with different CRSs (e.g., using st_distance()).
CRS Selection: Geographic CRS is best for global and spherical geometry operations. Projected CRS is ideal for local/regional planar geometry operations and analyses requiring precise distance and area calculations. Table 1 summarizes when to use geographic versus projected CRSs:
Code
library(tibble)library(gt)# Create the table as a tibbledata <-tibble(ID =1:9,Scenario =c("Global Mapping","Great Circle Distances","Angular Bearings","Radius-based Proximity Queries (Spherical)","Local Area Analysis","Buffer Creation","Overlay and Intersection Operations","Planar Geometry Operations","Area and Length Calculation" ),`Geographic CRS`=c("Required for displaying data on web maps (e.g., with leaflet).","Use for accurate calculation of shortest paths on the Earth's surface.","Essential for computing directions or angles between points.","Use for radius-based queries or geofencing on a global scale.","Not typically recommended unless working with global-scale spherical models.","Not recommended for smooth edges.","Limited to angular-based operations (e.g., great circle intersections).","Not ideal for creating planar geometries like buffers or centroids.","Use for global-scale spherical areas (with corrections)." ),`Projected CRS`=c("Not recommended due to distortions at global scales.","Distances may be inaccurate due to planar projections.","Distortions in angles make it unsuitable for such tasks.","Results can be distorted for large distances.","Ideal for small-to-medium-scale analyses with minimal distortions.","Best for creating accurate and smooth planar buffers.","Preferred for planar intersections and overlays for precise alignment.","Essential for tasks like buffering and centroids in local projections.","Use for precise calculations in local or regional contexts." ))# Create a beautiful table with gtdata |>gt() |>tab_header(title ="Geographic vs Projected CRS",subtitle ="A comparison of scenarios for applying Geographic and Projected CRSs" ) |> gtExtras::gt_theme_pff() |>tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns = Scenario) ) |>tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =`Geographic CRS`, # Column 3rows =1:4# Rows 1-4 ) ) |>tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =`Projected CRS`, # Column 4rows =5:9# Rows 5-9 ) )
Table 1: Guidelines for Choosing Between Geographic and Projected CRSs
Geographic vs Projected CRS
A comparison of scenarios for applying Geographic and Projected CRSs
ID
Scenario
Geographic CRS
Projected CRS
1
Global Mapping
Required for displaying data on web maps (e.g., with leaflet).
Not recommended due to distortions at global scales.
2
Great Circle Distances
Use for accurate calculation of shortest paths on the Earth's surface.
Distances may be inaccurate due to planar projections.
3
Angular Bearings
Essential for computing directions or angles between points.
Distortions in angles make it unsuitable for such tasks.
4
Radius-based Proximity Queries (Spherical)
Use for radius-based queries or geofencing on a global scale.
Results can be distorted for large distances.
5
Local Area Analysis
Not typically recommended unless working with global-scale spherical models.
Ideal for small-to-medium-scale analyses with minimal distortions.
6
Buffer Creation
Not recommended for smooth edges.
Best for creating accurate and smooth planar buffers.
7
Overlay and Intersection Operations
Limited to angular-based operations (e.g., great circle intersections).
Preferred for planar intersections and overlays for precise alignment.
8
Planar Geometry Operations
Not ideal for creating planar geometries like buffers or centroids.
Essential for tasks like buffering and centroids in local projections.
9
Area and Length Calculation
Use for global-scale spherical areas (with corrections).
Use for precise calculations in local or regional contexts.
7.6 Which CRS to use?
No single ‘best’ CRS: All projections involve distortions; the choice depends on the task.
Switch CRSs for different tasks: Use a CRS that aligns with your analysis or visualization goals.
Geographic CRS: most commonly, use WGS84 (EPSG:4326):
Common for web mapping, GPS datasets, and most vector - raster datasets.
EPSG code: 4326 (useful for converting to a universal CRS).
Projected CRS: Often dictated by local mapping agencies.
To ensure compatibility: Preferable to work in CRS (in which the data was provided)
UTM (Universal Transverse Mercator): a set of Projected CRSs
Note
The Universal Transverse Mercator (UTM) set of projections
Feature
Description
Projection Name
Universal Transverse Mercator (UTM)
Structure
Divides the Earth into 60 longitudinal zones (6° wide each) and 20 latitudinal segments.
CRS Type
Projected CRS. Uses the Transverse Mercator projection for each zone.
EPSG Codes
Northern Hemisphere: 32601 to 32660.
Southern Hemisphere: 32701 to 32760.
Zone Numbering
Zones numbered 1 to 60 (longitude).
Latitude segments labeled A to V (excluding I and O).
Example Zones
Zone 60H: Northern New Zealand, EPSG: 32760.
Zone 30N: Covers London, EPSG: 32630.
Key Features
Conformal projection: Preserves angles and shapes locally.
Suitable for small areas.
Distortion
Increases with distance from the zone’s central meridian.
Use Cases
Small-scale mapping, local surveys, engineering, and construction projects.
Reprojection: Transforming the coordinates of vector geometries (points, lines, polygons). Example: Converting CRS of an object for accurate distance measurement.
Key Functionality in sf Package: Use st_transform() to transform CRS.
Output contains components like User input and wkt.
CRS objects are lists containing elements like Name, proj4string, epsg, and $wkt, as shown in table in Section 7.3 above.
$wkt serves as the ultimate source of CRS details, queried via PROJ.
Note
OGR Coordinate Reference Systems and Coordinate Transformation tutorial provides insights into handling coordinate reference systems (CRS) and performing transformations using GDAL’s OGR library. Here’s how key concepts from the tutorial translate to {sf} in R:
Defining a Geographic CRS:
In {sf}, you can assign a CRS to an sf object using the st_set_crs() function.
Querying CRS Information:
In {sf}, you can retrieve CRS information using the st_crs() function:
Transforming Coordinates Between CRSs:
In {sf}, the st_transform() function is used to transform geometries to a different CRS. For example, to transform to EPSG:27700:
Handling Axis Order:
In {sf}, you can control axis order interpretation using the st_axis_order() function. To set it to traditional GIS order (longitude, latitude):
st_axis_order(FALSE)
7.8 Reprojecting Raster Geometries
Raster Reprojection: Creates a new raster with potentially different dimensions. Involves two steps:
Recommendation: Prefer vector reprojection over raster reprojection when both data types are used.
Key Difference from vector re-projection: Raster reprojection creates a new raster with modified attributes (e.g., columns, rows) rather than altering individual pixel coordinates like vectors.
Two Processes of raster reprojection:
Transformation: Adjusts grid geometry without changing pixel values, performed with the stars package. Also called Extent Transformation: Similar to vector reprojection, adjusts raster extent. Can be done using {stars} package. (Pebesma and Bivand 2023b)
Warping: Reprojects raster using methods like project() from the terra package. Also called Pixel Resampling: Computes new pixel values using methods like nearest neighbor or bilinear re-sampling. Can be done with terra::project() from {terra} (Hijmans 2024b)
Categorical Raster: Method used is nearest neighbor which ensures categorical values remain unchanged during reprojection. Result: Changes in resolution, extent, and added NA values (not new categories). Example is shown in Figure 3 .
Figure 3: Warping or Re-sampling categorical rasters (akin to reprojection) using terra::project() which uses near method - estimates re-sampled raster values using the value from the nearest original cell.
Continuous Raster: Method used is bilinear, computes pixel values as weighted averages of neighboring cells. Result: Modifications to mean values, resolution, and extent. Example is shown in Figure 4 .
Figure 4: Warping or Re-sampling a continuous value raster (akin to reprojection) using terra::project() which uses bilinear method - estimates re-sampled raster values using distance weighted mean from nearest 4 cells.
Projection Selection: Depends on the task. Use equal-area projections for density calculations, such as points per grid cell or inhabitants per grid cell.
Figure 5: Warping or Re-sampling elevation data for Greenland (obtained from {geodata} package) using terra::project() which uses bilinear method - estimates re-sampled raster values using distance weighted mean from nearest 4 cells.
7.9 Custom Map Projections
Established CRSs (e.g., EPSG:4326) are widely used but may not always fit specific needs; custom CRSs are an alternative.
Custom CRS Creation: Modify existing WKT CRS definitions, example shown in Figure 6, with:
For vector objects: sf::st_crs() and sf::st_transform()
For raster objects: terra::crs() and terra::project()
Code
# Example creating a custom azimuthal equidistant (AEQD) CRS projection focussed on Greenland# Deduce the centre of Greenlandgreenland |>st_centroid() |>st_geometry() |>st_as_text()# "POINT (-41.81367 73.08122)"st_crs(greenland)g1 <- greenland |>ggplot() +geom_sf() +labs(title ="WGS 84 (Geographic) CRS",subtitle ="Default Geographic CRS in the data" )custom_crs <- ('PROJCS["Custom_AEQD", GEOGCS["GCS_WGS_1984", DATUM["WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Azimuthal_Equidistant"], PARAMETER["Central_Meridian",-41.81367], PARAMETER["Latitude_Of_Origin",73.08122], UNIT["Meter",1.0]]')g2 <- greenland |>ggplot() +geom_sf() +geom_sf(data =st_centroid(greenland), size =1.5) +coord_sf(crs = custom_crs) +labs(title ="Custom AEQD CRS",subtitle ="Focussed around the centroid of Greenland (black dot)" )g <- g1 + g2 +plot_layout(nrow =1 ) +plot_annotation(title ="Creating a custom AEQD CRS for Greenland",tag_levels ="a",tag_prefix ="(",tag_suffix =")",theme =theme(plot.title =element_text(hjust =0.5, size =30) ) ) &theme(plot.title.position ="plot",plot.tag =element_text(size =24, face ="bold" ) )ggsave(plot = g,filename = here::here("book_solutions", "images", "chapter7-9_1.png"),height =700,width =1400,unit ="px",bg ="white")
Figure 6: An example creating a custom azimuthal equidistant (AEQD) CRS projection focussed on Greenland
Projection Wizard:
Interactive tool for selecting spatial extents and distortion properties (Projection Wizard).
Mollweide (+proj=moll): Preserves area relationships.
Winkel Tripel (+proj=wintri): Minimizes distortion across spatial properties.
Lambert Azimuthal Equal-Area (+proj=laea): Centered with +lon_0 and +lat_0 parameters. The parameters can be changed to suit the location you want to focus on, as shown in example below.
Code
# An example to plot the World map focussed around 4 major cities# Getting the World Mapworld_map <- rnaturalearth::ne_countries(scale ="small")# Create a tibble of major cities with their latitude and longitudemajor_cities <- tibble::tibble(city =c("New York", "London", "Tokyo", "Sydney"),latitude =c(40.7128, 51.5074, 35.6895, -33.8688),longitude =c(-74.0060, -0.1278, 139.6917, 151.2093))cities_geometry <- major_cities |>st_as_sf(coords =c("longitude", "latitude")) |>st_set_crs(value ="EPSG:4326")custom_plot <-function(city_name){# Select the city city_df <- major_cities |>filter(city == city_name)# Create custom LAEA projection focussed around the city coordinates new_crs =paste0("+proj=laea +x_0=0 +y_0=0 +lon_0=", city_df$longitude, " +lat_0=", city_df$latitude)ggplot() +geom_sf(data = world_map) +geom_sf(data = cities_geometry |>filter(city == city_name),size =2, colour ="red" ) +coord_sf(crs = new_crs) +labs(title = city_name,subtitle =paste0("Custom LAEA projection centered around ", city_name) ) +theme(axis.text =element_blank() )}g1 <-ggplot() +geom_sf(data = world_map) +labs(title ="World Map",subtitle ="Default Geographic CRS: WGS 84" )g6 <-ggplot() +geom_sf(data =st_transform(world_map, crs ="ESRI:54009")) +labs(title ="World Map (ESRI:54009)",subtitle ="Mollweide Projection" )g2 <-custom_plot("New York")g3 <-custom_plot("London")g4 <-custom_plot("Tokyo")g5 <-custom_plot("Sydney")my_design <- ("AAFFBCDE")g <-wrap_plots(g1, g2, g3, g4, g5, g6) +plot_layout(design = my_design,heights =c(1.4, 1) ) +plot_annotation(title ="Creating custom projections with {sf} in R",subtitle ="Lambert azimuthal equal-area projection, customized to be centered on any City around the globe",tag_levels ="a",tag_prefix ="(",tag_suffix =")",theme =theme(plot.title =element_text(hjust =0.5, size =60, face ="bold"),plot.subtitle =element_text(hjust =0.5, size =30, face ="bold") ) ) &theme(plot.title.position ="plot",plot.tag =element_text(size =24, face ="bold" ),plot.margin =margin(0,0,0,0, "pt") )ggsave(plot = g,filename = here::here("book_solutions", "images", "chapter7-9_2.png"),height =1400,width =2000,unit ="px",bg ="white")
Figure 7: An array of maps to demonstrate the power of R and {sf} to create custom projections. Part (a) shows the world map in Geographic CRS, best used for web-mapping. Part (b) shows the Mollweide Projection, showing the global map. The parts (c), (d), (e) and (f) are world maps focussed around a city around the globe in the Lambert Azimuthal Equal-Area projection, created using a custom function with coord_sf() and "+proj" strings.
Other Resources:
Learn about custom CRS definitions with WKT strings (OGC documentation).
Create a new object called nz_wgs by transforming nz object into the WGS84 CRS.
Create an object of class crs for both and use this to query their CRSs.
With reference to the bounding box of each object, what units does each CRS use?
Remove the CRS from nz_wgs and plot the result: what is wrong with this map of New Zealand and why?
The code chunk below shows the method to transform the CRS on nz object into WGS84 and querying their CRS using st_crs()$WktPretty |> str_view(). The New Zealand Transverse Mercator 2000 projection uses 1 metre as a unit, and WGS84 uses 1 degree of longitude / latitude as a unit. The result of removing the CRS from the nz_wgs object is plotted, along with original nz_wgs object, in Figure 8 .
g1 <-ggplot(data = nz_wgs) +geom_sf() +labs(title ="WGS84 Projection (Geographic)",subtitle ="The map is shown in aspect ratio as per the projection." ) +theme(plot.title.position ="plot" )nz_wgs1 <- nz_wgsst_crs(nz_wgs1) <-NAg2 <-ggplot(data = nz_wgs1) +geom_sf() +labs(title ="Same Map without projection (NA)",subtitle ="The map loses its fixed aspect ratio, and converts to equal axes spacing." ) +theme(plot.title.position ="plot" )g <- g1 + g2 +plot_layout(nrow =1 ) +plot_annotation(title ="Removing the CRS from New Zealand Map",tag_levels ="a",tag_prefix ="(",tag_suffix =")",theme =theme(plot.title =element_text(hjust =0.5, size =30) ) ) &theme(plot.title.position ="plot",plot.tag =element_text(size =24, face ="bold" ) )ggsave(plot = g,filename = here::here("book_solutions", "images", "chapter7-ex1.png"),height =700,width =1400,unit ="px",bg ="white")
Figure 8: Removing the CRS from a map causes coord_sf() to use coord_equal() instead of desirable aspect ratio prescribed in the original CRS displays.
E2.
Transform the world dataset to the transverse Mercator projection ("+proj=tmerc") and plot the result. What has changed and why? Try to transform it back into WGS 84 and plot the new object. Why does the new object differ from the original one?
The Figure 9 shows the entire result of transforming a Geographic CRS: map (a) in Figure 9, into Transverse Mercator Projection: map (b) in Figure 9, and then re-projecting it back into the same geographic CRS WGS84: map (c) in Figure 9 .
While geographic (GCS) and projected (PCS) systems allow versatile transformations, they are not perfectly reversible in real-world data due to floating-point precision limits and approximations inherent in projection algorithms. After transforming back to WGS 84, some distortions may appear because the process involves approximations during the transformation, especially if the intermediate PCS introduced distortion.
Code
world <- rnaturalearth::ne_countries(scale ="small", returnclass ="sf" ) |>select(name, geometry)# Removing Indonesia for easy plotting# filter(name != "Indonesia")g1 <-ggplot(data = world) +geom_sf() +labs(title ="Original World Map",subtitle ="Geographic CRS: WGS84" )g2 <-ggplot(data = world |>st_transform("+proj=tmerc")) +geom_sf() +labs(title ="Transverse Mercator Projection",subtitle ="Using st_crs(\"proj=tmerc\")" )g3 <-ggplot(data = world |>st_transform("+proj=tmerc") |>st_transform("WGS84") ) +geom_sf() +labs(title ="Back to original CRS (WGS84)",subtitle ="Russian and Antarctica segments are distorted." )g <- g1 + g2 + g3 +plot_layout(nrow =1 ) +plot_annotation(title ="Transforming, and reverting back from a Transverse Mercator to Geographic CRS",tag_levels ="a",tag_prefix ="(",tag_suffix =")",theme =theme(plot.title =element_text(hjust =0.5, size =42) ) ) &theme(plot.title.position ="plot",plot.tag =element_text(size =24, face ="bold" ) )ggsave(plot = g,filename = here::here("book_solutions", "images", "chapter7-ex2.png"),height =700,width =2400,unit ="px",bg ="white")
Figure 9: The entire result of transforming a Geographic CRS: original map (a), into (b) Transverse Mercator Projection, and then (c) back into the same geographic CRS WGS84.
E3.
Transform the continuous raster (con_raster) into NAD83 / UTM zone 12N using the nearest neighbor interpolation method. What has changed? How does it influence the results?
The code below, with the results shown in Figure 10, depicts an original continuous raster (con_raster) in part (a) of Figure 10, then the raster transformed into NAD83 / UTM Zone 12N with nearest neighbour method in part (b) of Figure 10, and finally, the raster transformed into NAD83 / UTM Zone 12N with bilinear method in part (c) of Figure 10 .
Changes:
The spatial alignment and resolution may change depending on the target CRS and grid size.
The edges of the raster may introduce distortions or artifacts, especially if the original CRS and UTM grid are significantly different.
Impact on Results:
Advantages: Nearest neighbor preserves the exact original values, which is critical for discrete data.
Disadvantages: For continuous rasters, it may introduce blocky artifacts, reducing the smoothness and potentially affecting analyses like slope or gradient computation.
Transform the categorical raster (cat_raster) into WGS 84 using the bilinear interpolation method. What has changed? How does it influence the results?
Transforming a categorical raster (cat_raster) from its original coordinate reference system (CRS) (NAD83 / UTM Zone 12N) (shown in (a) of Figure 11) into another CRS (such as WGS 84) using bilinear interpolation significantly alters its characteristics and interpretation, as shown in (b) of Figure 11.
Conversion of Raster Values:
A categorical raster represents discrete classes, such as land cover types (e.g., forest, urban, water), using integer values that act as category identifiers.
When transformed using bilinear interpolation, the algorithm treats these integers as continuous numeric values. It calculates a weighted average of the values from the nearest four cells in the input raster to assign values to the output raster.
This process results in non-integer values in the output raster, which no longer correspond to the original categories.
Loss of Categorical Meaning: The integer values in the original raster have a specific meaning tied to the categories they represent (e.g., 1 = forest, 2 = urban). After bilinear interpolation, the resulting non-integer values are meaningless in the context of categorical data because they no longer map directly to the predefined classes.
Blurring of Boundaries: Bilinear interpolation introduces smoothing between adjacent cells, which creates gradual transitions between values. This results in blurred boundaries between distinct categories, further distorting the original discrete nature of the data.
Introduction of Artefacts: The process can introduce artificial gradients or intermediate values that did not exist in the original data. For example, if the original raster has categories 1 (forest) and 2 (urban), interpolation might produce values like 1.5, which have no real-world significance.
This influences the results in the following ways: —
Inappropriate Representation: The output raster is now a continuous raster, which is unsuitable for applications requiring categorical data, such as land cover analysis or habitat classification. This could lead to misinterpretations if the data is used without understanding the transformation process.
Loss of Original Information: The categorical structure and associated metadata (e.g., legends, class descriptions) are effectively destroyed, rendering the raster unusable for its original purpose.
Impact on Analysis: Any analysis that relies on the classification of discrete land cover types (e.g., area statistics, habitat mapping) would yield invalid results, as the categories no longer exist in their original form.
Thus, for transforming categorical rasters between CRSs, nearest-neighbor interpolation (show in part (c) of Figure 11) should be used instead of bilinear interpolation. Nearest-neighbor interpolation preserves the integer values of the categories, ensuring that the categorical meaning of the data is retained.
Code
cat_raster <-rast(system.file("raster/nlcd.tif", package ="spDataLarge"))# Get the CRS of the rastercrs(cat_raster) |>str_view()# Check associated colours with the raster# has.colors(cat_raster)# coltab(cat_raster)# My custom palette for uniformitymypal <-c("#008ee0", "darkred", "grey80", "darkgreen","#858f42", "#90bf7e", "#d9e86b", "#96defa")g1 <-ggplot() +geom_spatraster(data = cat_raster) +scale_fill_manual(values = mypal) +coord_sf(expand =FALSE) +labs(title ="Categorical Raster",subtitle ="CRS NAD83 / UTM Zone 12 N",fill ="Land Cover Type" )g2 <-ggplot() +geom_spatraster(data = cat_raster |>project("WGS84", method ="bilinear") ) +coord_sf(expand =FALSE) +scale_fill_princess_c(na.value ="black") +labs(title ="Projected Raster (bilinear)",subtitle ="Raster converted into a Continuous Raster.\nCategories are Lost !",fill ="Value" )g3 <-ggplot() +geom_spatraster(data = cat_raster |>project("WGS84", method ="near") ) +coord_sf(expand =FALSE) +scale_fill_manual(values = mypal, na.value ="black") +labs(title ="Projected Raster (nearest neighbour)",subtitle ="Interpolation into Geographic CRS: WGS84.\nA new category of `NA` is introduced, as expected.",fill ="Land Cover Type" )g <- g1 + g2 + g3 +plot_layout(nrow =1 ) +plot_annotation(title ="Transforming a categorical raster using the nearest neighbor and bilinear interpolation method.",subtitle ="Bilinear interpolation computes averages, and leads the raster to lose all meaning, and convert into a continuous raster. (NA values are in black).",tag_levels ="a",tag_prefix ="(",tag_suffix =")",theme =theme(plot.title =element_text(hjust =0.5, size =36),plot.subtitle =element_text(hjust =0.5, size =24) ) ) &theme(plot.title.position ="plot",plot.tag =element_text(size =24, face ="bold" ),legend.key.width =unit(5, "pt"),plot.margin =margin(0,5,0,5, "pt"),legend.margin =margin(0,0,0,0, "pt"),legend.box.margin =margin(0,0,0,0, "pt") )ggsave(plot = g,filename = here::here("book_solutions", "images", "chapter7-ex4.png"),height =800,width =2100,unit ="px",bg ="white")
Figure 11
References
Cheng, Joe, Barret Schloerke, Bhaskar Karambelkar, and Yihui Xie. 2024. “Leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library.”https://CRAN.R-project.org/package=leaflet.
