0Pricing
R Academy · Lesson

Coordinate Reference Systems and Projections

Transform between CRS using st_crs() and st_transform() correctly.

Coordinate Reference Systems and Projections is a free R Academy lesson on CoddyKit — lesson 2 of 4. You can read the complete lesson below for free — then practise it hands-on in the browser with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, one of 4 lessons in the course, and your progress syncs across the web and the CoddyKit app.

Why CRS Matters

A Coordinate Reference System (CRS) tells R what the numbers in your geometry column mean — degrees of latitude/longitude, metres in a projected grid, or something else. Mixing data in different CRS without re-projecting causes spatial joins to fail silently and distance calculations to be wrong.

library(sf)

# Two sf objects in the same CRS: operations work
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))

# sf will error or warn when CRS do not match during operations
# This is the CRS metadata attached to the object
crs <- st_crs(nc)
cat('EPSG:', crs$epsg, '\n')
cat('Input is geographic (lon/lat)?', st_is_longlat(nc), '\n')
cat('Units:', crs$units_gdal, '\n')

st_crs: Reading CRS Info

st_crs(sf_obj) returns a crs object with EPSG code, WKT string, Proj4 string, and unit information. You can also pass an integer EPSG code directly to retrieve the CRS definition.

library(sf)

nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))

# Inspect the CRS object
crs <- st_crs(nc)
cat('Class:', class(crs), '\n')
cat('EPSG:', crs$epsg, '\n')
cat('Proj4:', substr(crs$proj4string, 1, 60), '...\n')

# Lookup a CRS by EPSG code
wgs84 <- st_crs(4326)
web_mercator <- st_crs(3857)

cat('\nWGS84 name:', wgs84$Name, '\n')
cat('Web Mercator name:', web_mercator$Name, '\n')

EPSG:4326 — WGS84

EPSG:4326 is the World Geodetic System 1984 — the standard for GPS coordinates. Coordinates are in decimal degrees: longitude (X) ranges −180 to 180, latitude (Y) ranges −90 to 90. It is a geographic CRS, not projected.

library(sf)

# Create points in WGS84 (GPS coordinates)
cities <- data.frame(
  city = c('London', 'Tokyo', 'Sydney', 'New York'),
  lon  = c(-0.118, 139.692, 151.209, -74.006),
  lat  = c(51.509,  35.689, -33.869,  40.713)
)

cities_sf <- st_as_sf(cities,
  coords = c('lon', 'lat'),
  crs    = 4326
)

cat('CRS:', st_crs(cities_sf)$Name, '\n')
cat('Units:', st_crs(cities_sf)$units_gdal, '\n')

# Naive distance in degrees (not meaningful for real distances)
d <- st_distance(cities_sf[1, ], cities_sf[2, ])
cat('London-Tokyo distance (degrees):', round(as.numeric(d), 2), '\n')

EPSG:3857 — Web Mercator

EPSG:3857 is the Web Mercator projection used by Google Maps, OpenStreetMap, and Leaflet. Coordinates are in metres, centred at the prime meridian / equator. It distorts areas near the poles but is good for web visualisation.

library(sf)

cities <- data.frame(
  city = c('London', 'Tokyo'),
  lon  = c(-0.118, 139.692),
  lat  = c(51.509,  35.689)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)

# Transform to Web Mercator
cities_merc <- st_transform(cities_sf, 3857)

cat('WGS84 coords (degrees):\n')
print(st_coordinates(cities_sf))

cat('\nWeb Mercator coords (metres):\n')
print(st_coordinates(cities_merc))

# Distance is now in metres
d <- st_distance(cities_merc[1, ], cities_merc[2, ])
cat('\nLondon-Tokyo distance (km):', round(as.numeric(d) / 1000, 0), '\n')

st_set_crs: Assigning a CRS

st_set_crs(sf_obj, crs) assigns a CRS without transforming coordinates. Use this only when the CRS metadata is missing or wrong. If you want to re-project coordinates to a new CRS, use st_transform() instead.

library(sf)

# Create geometry with NO CRS
pts <- st_sfc(
  st_point(c(-74.006, 40.713)),
  st_point(c(-87.629, 41.878))
)
cat('CRS before:', is.na(st_crs(pts)), '(NA = no CRS)\n')

# Assign WGS84 (no coordinate transformation)
pts_wgs84 <- st_set_crs(pts, 4326)
cat('CRS after assignment:', st_crs(pts_wgs84)$epsg, '\n')

# Compare with st_transform (DOES transform coordinates)
pts_sf    <- st_sf(id = 1:2, geometry = pts_wgs84)
pts_merc  <- st_transform(pts_sf, 3857)
cat('Original lon:', st_coordinates(pts_sf)[1, 'X'], '\n')
cat('Mercator X:  ', st_coordinates(pts_merc)[1, 'X'], '\n')

st_transform: Re-Projecting

st_transform(sf_obj, crs) re-projects all geometries from their current CRS to the target CRS, applying the full geodetic transformation. The crs argument accepts an EPSG integer, a Proj4 string, or a WKT string.

library(sf)

nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
cat('Original CRS:', st_crs(nc)$epsg, '\n')

# Re-project to WGS84
nc_wgs84 <- st_transform(nc, 4326)
cat('Transformed CRS:', st_crs(nc_wgs84)$epsg, '\n')

# Re-project to a local projected CRS (NC State Plane)
nc_stateplane <- st_transform(nc, 32119)  # NAD83 / North Carolina
cat('State Plane CRS:', st_crs(nc_stateplane)$Name, '\n')

# Area calculation is more accurate in a projected CRS
area_deg <- st_area(nc[1, ])
area_m2  <- st_area(nc_stateplane[1, ])
cat('Area (degrees^2):', format(area_deg, big.mark = ','), '\n')
cat('Area (m^2):      ', format(area_m2,  big.mark = ','), '\n')

Choosing the Right Projection

Different projections optimise for different properties: equal-area (e.g. Mollweide), conformal (Mercator), equidistant. For local analysis use a projected CRS in metres; for global overviews WGS84 works; for web maps use Web Mercator.

library(sf)

# Common EPSG codes cheat sheet
epsg_table <- data.frame(
  EPSG = c(4326, 3857, 4269, 32637, 27700),
  Name = c(
    'WGS84 (GPS, global)',
    'Web Mercator (web maps)',
    'NAD83 (North America geographic)',
    'UTM zone 37N (Middle East)',
    'British National Grid (UK)'
  ),
  Units = c('degrees', 'metres', 'degrees', 'metres', 'metres')
)

print(epsg_table, row.names = FALSE)

# Look up any EPSG online: https://epsg.io/<code>
cat('\nTip: use st_crs(epsg_code)$Name to verify a code\n')

Proj4 Strings

Before EPSG codes became standard, CRS were defined by Proj4 strings — compact parameter strings passed to the PROJ library. You can still use them with st_transform(), but EPSG codes are preferred for clarity and future compatibility.

library(sf)

nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))

# WGS84 as Proj4 string
wgs84_proj4 <- '+proj=longlat +datum=WGS84 +no_defs'
nc_wgs84 <- st_transform(nc, crs = wgs84_proj4)

cat('Using Proj4 string:\n', wgs84_proj4, '\n')
cat('Result CRS EPSG:', st_crs(nc_wgs84)$epsg, '\n')

# Compare: equivalent to st_transform(nc, 4326)
cat('Same as EPSG 4326?',
    identical(st_crs(nc_wgs84), st_crs(st_transform(nc, 4326))), '\n')

Checking CRS Equality

Before any spatial join or overlay operation, verify that both objects share the same CRS using st_crs(a) == st_crs(b). If they differ, transform one to match the other before proceeding.

library(sf)

# Two objects in different CRS
nc_nad27  <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84  <- st_transform(nc_nad27, 4326)

cat('Same CRS?', st_crs(nc_nad27) == st_crs(nc_wgs84), '\n')

# Always check before join/overlay
align_crs <- function(a, b) {
  if (st_crs(a) != st_crs(b)) {
    message('CRS mismatch! Transforming b to match a.')
    b <- st_transform(b, st_crs(a))
  }
  list(a = a, b = b)
}

aligned <- align_crs(nc_nad27, nc_wgs84)
cat('After alignment - same CRS?',
    st_crs(aligned$a) == st_crs(aligned$b), '\n')

Calculating Distances with CRS

Distances computed on geographic CRS (degrees) are meaningless. Re-project to an equal-distance or equal-area projected CRS first, or use sf's built-in geodesic distance calculation which handles WGS84 correctly via the PROJ library.

library(sf)

# City points in WGS84
cities <- data.frame(
  city = c('Paris', 'Berlin', 'Madrid', 'Rome'),
  lon  = c(2.350, 13.404, -3.702, 12.496),
  lat  = c(48.865, 52.520,  40.417, 41.902)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)

# sf computes geodesic distances on WGS84 automatically
dist_matrix <- st_distance(cities_sf)
rownames(dist_matrix) <- cities$city
colnames(dist_matrix) <- cities$city

# Convert to km
dist_km <- round(as.numeric(dist_matrix) / 1000)
mat_km <- matrix(dist_km, 4, 4, dimnames = list(cities$city, cities$city))
print(mat_km)

A Complete CRS Workflow

Here is a complete workflow: load raw GPS data, assign WGS84, re-project to a local metric CRS for accurate area/distance calculations, perform the analysis, and re-project back to WGS84 for web mapping output.

library(sf)

# Load and assign CRS
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84 <- st_transform(nc, 4326)

# Re-project to NC State Plane (metres) for accurate areas
nc_proj <- st_transform(nc, 32119)

# Compute area per county in km2
nc_proj$area_km2 <- as.numeric(st_area(nc_proj)) / 1e6

# Top 5 largest counties
top5 <- nc_proj[order(-nc_proj$area_km2), c('NAME', 'area_km2')]
cat('5 Largest NC Counties (km2):\n')
print(head(as.data.frame(top5)[, c('NAME', 'area_km2')], 5))

# Back to WGS84 for output
nc_output <- st_transform(nc_proj, 4326)
cat('\nOutput CRS:', st_crs(nc_output)$epsg, '\n')

Quick Check

You need to compute accurate areas (in square metres) for polygon features currently in WGS84 (EPSG:4326). What should you do?

Recap: CRS and Projections

Key takeaways:

  • st_crs(sf_obj) reads CRS metadata; EPSG codes are the preferred way to specify CRS
  • EPSG:4326 = WGS84 (degrees, for GPS); EPSG:3857 = Web Mercator (metres, for web maps)
  • st_set_crs() assigns CRS without transforming; st_transform() re-projects coordinates
  • Always verify CRS equality with st_crs(a) == st_crs(b) before spatial operations
  • For accurate distances and areas, re-project to a metric projected CRS first
  • Proj4 strings still work but EPSG codes are cleaner and more readable
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))

# CRS pipeline
cat('Original:', st_crs(nc)$epsg, '\n')
nc_wgs84 <- st_transform(nc, 4326)
nc_proj  <- st_transform(nc_wgs84, 32119)
cat('WGS84:', st_crs(nc_wgs84)$epsg, '\n')
cat('Projected:', st_crs(nc_proj)$epsg, '\n')

Frequently asked questions

Is the “Coordinate Reference Systems and Projections” lesson free?

Yes — the full text of “Coordinate Reference Systems and Projections” is free to read here on the web, and the R Academy course includes 4 lessons in total. To practise it interactively (a built-in code editor and a 24/7 AI tutor) and unlock the rest of the R Academy course, upgrade to CoddyKit PRO.

What will I learn in “Coordinate Reference Systems and Projections”?

Transform between CRS using st_crs() and st_transform() correctly. You practise R Academy with hands-on code you run directly in the browser, and a 24/7 AI tutor answers your questions as you work through the lesson.

Do I need any experience to start R Academy?

No prior experience is required. R Academy on CoddyKit is structured for beginners through advanced learners; this is — lesson 2 of 4, so you can start here or from the beginning and move at your own pace.

How long does the “Coordinate Reference Systems and Projections” lesson take?

Most CoddyKit lessons take about 5–10 minutes. Each one is bite-sized and interactive, so you make steady progress and pick up exactly where you left off across the web and the app.

Can I write and run code in this R Academy lesson?

Yes. Every R Academy lesson includes a built-in code editor, so you write and run real code right in your browser and get instant AI feedback — no local setup required.

All lessons in this course

  1. Simple Features with the sf Package
  2. Coordinate Reference Systems and Projections
  3. Spatial Joins and Operations
  4. Interactive Maps with leaflet
← Back to R Academy