axis order in GDAL >= 2.5.0

[edzer]

Should sf worry about axis order? See rfc 73.

[rsbivand]

I think that sp-facing code stays with GIS-order (easting, northing), and I'll have to fix rgdal to stabilise that. For sf, the established revdep codebase is smaller, and we know much less about future incoming data. If future incoming data will expect axis order to be respected, we may need to handle that, converting to GIS-order late (geom_sf(), tmap::tm_shape(), mapview::mapview()) based on the CRS, which would be the valid source of authority. We'd need PROJ/GDAL to provide tools accessing the axis order for a given CRS to know whether to late-swap. Are those tools there (SetAxisMappingStrategy())? Do we want stars/sf to only write GIS-order files? If proxy data sources are not GIS-order, late-swapping for visualisation may make sense. Could we survey people - what is leaflet going to do?

[edzer]

I'm still trying to get my head around the simplest things, like why with gdal 2.5.0dev1

> st_transform(st_sfc(st_point(0:1), crs = 4326), 3857)
Geometry set for 1 feature 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0 ymin: 111325.1 xmax: 0 ymax: 111325.1
epsg (SRID):    3857
proj4string:    +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
POINT (0 111325.1)

did not assume lat/long, and

> st_crs(4326)
Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

[rsbivand]

Indeed, for 4326 EPSG itself says:

Axes: latitude, longitude.

on http://www.epsg-registry.org/ querying code 4326.

[edzer]

Example in r-spatial/discuss#28

[edzer]

Hoho! We can now swap axes:

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
st_sfc(p1, p2, crs = 4326) %>% 
  st_transform(pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1")
# Geometry set for 2 features 
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 20 ymin: -30 xmax: 52 ymax: 7
# geographic CRS: WGS 84
# epsg (SRID):    4326
# POINT (52 7)
# POINT (20 -30)

[edzer]

This I now can get my head around:

library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
# load the nc dataset, EPSG 4326, which is stored with long-lat coordinates:
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
st_bbox(nc)
#      xmin      ymin      xmax      ymax 
# -84.32385  33.88199 -75.45698  36.58965 
# swap the coordinates to lat-long, respecting EPSG 4326 authority:
nc = st_transform(nc, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1")
# verify:
st_bbox(nc) # x=first, y=second
#      xmin      ymin      xmax      ymax 
#  33.88199 -84.32385  36.58965 -75.45698 
# instrument sf such that it respects EPSG authority axis order (i.e., lat long for 4326):
st_axis_order(authority_compliant = TRUE)
# now transform to NC state plane:
nc.sp = st_transform(nc, 2264)
# plot:
plot(nc.sp[,1], axes = TRUE, graticule = TRUE) # good

with the commit to plot.sf, axes are swapped back on the fly when needed for plotting. (This seems much more tricky for the plot.sfc_* methods.)

I guess this can be closed now.

[rsbivand]

I did speculate in https://r-spatial.github.io/spdep/articles/sids.html that sids is actually +proj=longlat +ellps=clrk66 +datum=NAD27; EPSG:4267:

> cat(showSRID("EPSG:4267", multiline="YES"), "\n")
GEOGCRS["NAD27",
    DATUM["North American Datum 1927",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            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["North America - NAD27"],
        BBOX[7.15,167.65,83.17,-47.74]],
    ID["EPSG",4267]] 

[edzer]

That is what I see in nc too:

> st_crs(nc)
Coordinate Reference System:
  User input: NAD27 
  wkt:
GEOGCRS["NAD27",
    DATUM["North American Datum 1927",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            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["North America - NAD27"],
        BBOX[7.15,167.65,83.17,-47.74]],
    ID["EPSG",4267]]