Skip to main content

A Minimal Set of Lectures Notes

6 Coordinate Reference Systems

6.1 References

6.1.1 Read Chapter 09 – Coordinate Systems

6.1.2 Chapter - Coordinate Reference Systems

6.1.3 Very brief recap of cartography:

6.1.3.2 Concepts

  • Meridians (Greenwich) and Parallels (Equator);
  • Geoid vs Ellipsoid (WGS84);
  • Coordinate Reference System or CRS (Geocentric, Geografic, Projection);
  • Trasformations from one CRS definition to another;
  • Datum (WGS84, Roma40, …):
  • Spatial statistical analysis requires projections!

6.3 Coordinates Conversion

6.3.1 Search in R help | transform and project 

??transform  # search for string 'transform' in R Help

We get several options. The following functions are selected because they are related to both vector

  • PROJ::proj_trans()
  • proj4::ptransform
  • sf::st_tranform()

and raster geodata:

  • stars::st_transform
??project # search for string 'project' in R Help

We get above all one option, which is more related to raster geodata:

  • raster::projectRaster()
  • terra::project()

6.3.2 Use R

The suggestion is to build the required data model (i.e. vector or raster) and then proceed with the transform function in R.

6.3.2.1 vector data model

Build a simple feature object, and then use st_transform():

# The point below is taken from Google Maps, which has EPSG:4326
x <- 14.516343  # Longitude
y <- 41.224083  # Latitude
google_crs <- 4326
require(sf)
#> Loading required package: sf
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2()
#> is TRUE
P_goo <- st_point(c(x, y))
P_goo
#> POINT (14.51634 41.22408)
P_goo_sfc <- st_sfc(P_goo, crs = google_crs)
class(P_goo_sfc)
#> [1] "sfc_POINT" "sfc"
P_goo_sfc
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 14.51634 ymin: 41.22408 xmax: 14.51634 ymax: 41.22408
#> Geodetic CRS:  WGS 84
#> POINT (14.51634 41.22408)
my_crs <- 32633
P_my_sf <- st_transform(P_goo_sfc, my_crs)
P_my_sf
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 459461.7 ymin: 4563746 xmax: 459461.7 ymax: 4563746
#> Projected CRS: WGS 84 / UTM zone 33N
#> POINT (459461.7 4563746)
P_goo_sf <- st_sf(point = 1, geometry = P_goo_sfc)
class(P_goo_sf)
#> [1] "sf"         "data.frame"
P_goo_sf$geometry
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 14.51634 ymin: 41.22408 xmax: 14.51634 ymax: 41.22408
#> Geodetic CRS:  WGS 84
#> POINT (14.51634 41.22408)
my_crs <- 32633
P_my_sf <- st_transform(P_goo_sf, my_crs)
P_my_sf$geometry
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 459461.7 ymin: 4563746 xmax: 459461.7 ymax: 4563746
#> Projected CRS: WGS 84 / UTM zone 33N
#> POINT (459461.7 4563746)

6.3.2.2 raster data model

Define the coordinates of the Upper-Left and Lower-Right points delimiting the bounding box of the RASTER we want to create:

# Click on Google maps on the two points of the bbox and write the coordinates values here
UL <- c(11, 43)  # upper left
LR <- c(17, 38)  # lower right
require(terra)
#> Loading required package: terra
#> terra 1.8.54
GRID <- rast( xmin = UL[1], 
              xmax = LR[1], 
              ymin = LR[2], 
              ymax = UL[2], 
              resolution = 0.1, 
              crs = "epsg:4326" )
GRID
#> class       : SpatRaster 
#> size        : 50, 60, 1  (nrow, ncol, nlyr)
#> resolution  : 0.1, 0.1  (x, y)
#> extent      : 11, 17, 38, 43  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
values(GRID) <- 1:ncell(GRID)
plot(GRID)
GRID_pr <- project(GRID, "EPSG:32633")

Plot: note that it is slightly rotated with respect to the previous one due to reprojection

plot(GRID_pr)

6.3.3 Use epsg.io web application

https://epsg.io/transform.

For instance, convert from 3226 into 3857.

s_srs <- 4326
t_srs <- 32633
x <- 14.516343  # Longitude
y <- 41.224083  # Latitude
example_URL <- "https://epsg.io/transform#s_srs=4326&t_srs=32633&x=13.5000000&y=42.5000000"
my_URL <- paste0("https://epsg.io/transform#s_srs=", s_srs, "&t_srs=", t_srs, "&x=", x, "&y=", y)
my_URL
#> [1] "https://epsg.io/transform#s_srs=4326&t_srs=32633&x=14.516343&y=41.224083"
xc <- 459461.72
yc <- 4563745.65

Vienna, Austria
P1(15.549644, 48.374276)
P2(17.527792, 47.615370)

P1 <- st_point(c(15.549644, 48.374276))
P2 <- st_point(c(17.527792, 47.61537))
P <- st_sfc(P1, P2, crs = 4326)
P_3035 <- st_transform(P, crs = 3035)
P_3035
#> Geometry set for 2 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 4731873 ymin: 2751202 xmax: 4886176 ymax: 2822213
#> Projected CRS: ETRS89-extended / LAEA Europe
#> POINT (4731873 2822213)
#> POINT (4886176 2751202)