Site icon R-bloggers

R spatial follows GDAL and PROJ development

[This article was first published on r-spatial, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

[view raw Rmd]

GDAL and PROJ

GDAL and PROJ (formerly proj.4) are two libraries that form the basis, if not foundations, for most open source geospatial software, including most R packages (sf, sp, rgdal, and all their dependencies). The dependency for package sf is for instance pictured here:

Briefly:

gdalbarn

Motivated by the need for higher precision handling of coordinate transformations and the wish to support for a better description of coordinate reference systems (WKT2), a succesful fundraising helped the implementation of a large number of changes in GDAL and PROJ, most notably:

crs objects in sf

Pre-0.9 versions of sf used crs (coordinate reference system) objects represented as lists with two components, epsg (possibly set as NA) and proj4string:

library(sf) 
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
st_crs(4326)
# Coordinate Reference System:
#   EPSG: 4326
#   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

now, with sf >= 0.9, crs objects are lists with two components, input and wkt:

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1

(x = st_crs(4326))

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

where a $ method allows for retrieving the epsg and proj4string values:

x$epsg

## [1] 4326

x$proj4string

## [1] "+proj=longlat +datum=WGS84 +no_defs"

but this means that packages that hard-code for instance

x[["proj4string"]]

## NULL

now fail to get the result wanted; NULL is not a value that would have occurred in legacy code.

Regretably, assignment to a crs object component still works, as the objects are lists, so not all downstream legacy code will fail

x$proj4string <- "+proj=longlat +ellps=intl"
x$proj4string

## Warning in `$.crs`(x, proj4string): old-style crs object found: please update
## code

## [1] "+proj=longlat +ellps=intl +no_defs"

Package maintainers and authors of production scripts will need to review their use of crs objects.

Many external data sources provide a WKT CRS directly and as such do not have an “input” field. In such cases, the input field is filled with the CRS name, which is a user-readable representation

st = stars::read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
st_crs(st)$input

## [1] "UTM Zone 25, Southern Hemisphere"

but this representation can not be used as input to a CRS:

st_crs(st_crs(st)$input)

## Error in st_crs.character(st_crs(st)$input): invalid crs: UTM Zone 25, Southern Hemisphere

however wkt fields obviously can be used as input:

st_crs(st_crs(st)$wkt) == st_crs(st)

## [1] TRUE

CRS objects in sp

When equiped with a new ( >= 1.5.6) rgdal version, sp’s CRS objects carry a comment field with the WKT representation of a CRS:

# install.packages("rgdal", repos="http://R-Forge.R-project.org")
library(sp)
(x = CRS("+init=epsg:4326")) # or better: CRS(SRS_string='EPSG:4326')

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

cat(comment(x), "\n")

## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]]]

and it is this WKT representation that is used to communicate with GDAL and PROJ when using packages rgdal or sf. At present, rgdal generates many warnings about discarded PROJ string keys, intended to alert package maintainers and script authors to the need to review code. It is particularly egregious to assign to the CRS object projargs slot directly, and this is unfortunately seem in much code in packages.

Coercion from CRS objects to crs and back

Because workflows often need to combine packages using sp and sf representations, coercion methods from CRS to crs have been updated to use the WKT information; from sp to sf one can use

(x2 <- st_crs(x))

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]]]

The sp CRS constructor has been provided with an additional argument SRS_string= which accepts WKT, among other representations

(x3 <- CRS(SRS_string = x2$wkt))

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

but also

(x4 <- as(x2, "CRS"))

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

uses the WKT information when present.

all.equal(x, x3)

## [1] TRUE

all.equal(x, x4)

## [1] TRUE

Axis order

R-spatial packages have, for the past 25 years, pretty much assumed that two-dimensional data are XY-ordered, or longitude-latitude. Geodesists, on the other hand, typically use \((\phi,\lambda)\), or latitude-longitude, as coordinate pairs; the PROJ logo is now PR\(\phi\)J. If we use geocentric coordinates, there is no logical ordering. Axis direction may also vary; the y-axis index of images typically increases when going south. As pointed out in sf/#1033, there are powers out there that will bring us spatial data with (latitude,longitude) as (X,Y) coordinates. Even stronger, officially, EPSG:4326 has axis order latitude, longitude (see WKT description above).

Package sf by default uses a switch in GDAL that brings everything in the old, longitude-latitude order, but data may come in in another ordering.

This can now be controlled (to some extent), as st_axis_order can be used to query, and set whether axis ordering is “GIS style” (longitude,latitude; non-authority compliant) or “authority compliant” (often: latitude,longitude):

pt = st_sfc(st_point(c(0, 60)), crs = 4326)
st_axis_order() # query default: FALSE means interpret pt as (longitude latitude)

## [1] FALSE

st_transform(pt, 3857)[[1]]

## POINT (0 8399738)

(old_value = st_axis_order(TRUE)) 

## [1] FALSE

# now interpret pt as (latitude longitude), as EPSG:4326 prescribes:
st_axis_order() # query current value

## [1] TRUE

st_transform(pt, 3857)[[1]]

## POINT (6679169 0)

st_axis_order(old_value) # set back to old value

sf::plot is sensitive to this and will swap axis if needed, but for instance ggplot2::geom_sf is not yet aware of this.

Workflows using sp/rgdal should expect “GIS style” axis order to be preserved

rgdal::get_enforce_xy()

## [1] TRUE

pt_sp <- as(pt, "Spatial")
coordinates(pt_sp)

##      coords.x1 coords.x2
## [1,]         0        60

coordinates(spTransform(pt_sp, CRS(SRS_string="EPSG:3857")))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +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

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

##      coords.x1 coords.x2
## [1,]         0   8399738

Further reading

To leave a comment for the author, please follow the link and comment on their blog: r-spatial.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.