WKT axis order and CRS:84

[dblodgett-usgs]

I hesitate to bring this up, but there's been some discussion of what the "right" thing to do is regarding axis order and use of EPSG:4326 with WKT geometry in a project I am taking part in.

Technically, when writing out WKT, if 4326 is included as the SRID, the axis order should reflect lat - lon. Practically, (and this is where I might get some 🔥) we should always use lon - lat axis order to avoid the confusion and misinterpretation introduced by the few clients that us the lat - lon order. I'm not going to suggest changing to y/x axis order for particular CRSs by default, but it would be nice to:

1. have optional "strict" treatment of axis order through a T/F flag on st_as_text() and
2. allow specification of a lon/lat CRS in the WKT string.

I'm not sure how this would work with the SRID=#### approach used with the WKT in sf as there isn't an EPSG code for the urn:ogc:def:crs:OGC:1.3:CRS84 but it would be really nice if there were a way to emit the WKT string with an explicitely lon/lat axis order CRS. Maybe the approach taken in GeoSPARQL described here would be an option?

so... some code to illustrate

This all works now:

lat_lon <- sf::st_sf(id = 1, geom = sf::st_sfc(sf::st_point(c(-90,40))))
sf::st_crs(lat_lon) <- sf::st_crs(4326)
sf::st_as_text(lat_lon$geom)
# [1] "POINT (-90 40)"
sf::st_as_text(lat_lon$geom, EWKT = TRUE)
# [1] "SRID=4326;POINT (-90 40)"

This is not the current behavior, but illustrates what might be done:

sf::st_as_text(lat_lon$geom, EWKT = TRUE, STRICT_AXIS_ORDER = TRUE)
# [1] "SRID=4326;POINT (40,-90)"
sf::st_as_text(lat_lon$geom, EWKT = TRUE, STRICT_AXIS_ORDER = FALSE)
# [1] "<http://www.opengis.net/def/crs/OGC/1.3/CRS84> POINT (-90, 40)"

That all said, I would be just as happy with the statement that the axis order of lat/lon CRSs is always lon - lat regardless of the SRID of the data. But a statement to that affect should probably be included in the doc for st_as_text()?

[edzer]

Major can of worms opening. Could you please elaborate a bit on three things:

1. which practical problem for you would be solved by this
2. where (or whether), in the SFA standard, the axis order of coordinate pairs is determined, and allowed to be lat long?
3. which spatial databases do this, and how?

[dblodgett-usgs]

Agreed -- can of worms. I only bring this up because this issue will not go away and I think it's in our interest to be up front about our handling.

1. First, and practically speaking, according to EPSG:4326, coordinate pairs should be in lat/lon (y/x) order. So, if I had to interpret the string: SRID=4326;POINT (-90 40) seeing SRID=4326 I would have to ask myself, "is this really lat/lon? Or is this the more typical (and not strictly conformant to 4326) lon/lat?" Second, for me personally, it would be helpful to have a WKT output that conforms to the WKT convention used by GeoSPARQL e.g. <http://www.opengis.net/def/crs/OGC/1.3/CRS84> POINT (-90, 40) or <http://www.opengis.net/def/crs/EPSG/0/4326> POINT (40, -90) By my understanding of the goals of sf this is likely not in the intended scope?
2. I've just spent a bit of time with the SFA spec and googled around for more. Indeed, SFA prescribes lon/lat axis order. This reference is very helpful and speaks to the WKT/EWKT conventions use of lon/lat order exclusively. Excerpting the pertinent line:

In the [SIMPLE-FEATURES] definition of WKT, the coordinate axis-order is by default longitude / latitude, irrespective of the coordinate reference system used. The same applies to EWKT (Extended WKT) - a PostGIS extension to WKT supported also by other GIS tools - which includes a parameter (SRID) for specifying the coordinate reference system. For this reason, whenever using WKT to encode geometries, it is important that the reference WKT specification can be unambiguously determined.

3. The best treatment of this issue that I'm aware of is in GeoTools. This is a great discussion of the issue.

At the end of the day, I think the best thing to do here is just declare what sf does. That is, add a line to the documentation of sf::st_as_text() stating that the WKT emitted conforms to the SFA WKT and the EWKT extension introduced by PostGIS and broadly implemented to support declaring SRIDs. I would be happy to prepare and submit a PR to this affect if it is agreeable.

[edzer]

Thanks for the clarifications and links to docs on this issue. I'd welcome your proposed PR!

[edzer]

#1033