importEPSG and importEPSGA both swap axis order to lat/lon

[rafaqz]

ArchGDAL.jl/src/spatialref.jl

Line 145 in f57c2dc

This method is similar to importFromEPSGA() except that EPSG preferred axis

I think this no longer happens as of GDAL 3.0, and both methods are identical:

https://gdal.org/doxygen/classOGRSpatialReference.html#aaa6965a1df98cdc673dfb20697eab613

julia> p = ArchGDAL.createpoint(0.0, 180.0)
Geometry: POINT (0 180)

julia> ArchGDAL.importEPSG(4326) do sourcecrs_ref
           ArchGDAL.importPROJ4("+proj=cea") do targetcrs_ref
               ArchGDAL.createcoordtrans(sourcecrs_ref, targetcrs_ref) do transform
                   ArchGDAL.transform!(p, transform)
               end
           end
       end
Geometry: POINT (20037508.3427892 0.0)

julia> p = ArchGDAL.createpoint(0.0, 180.0)
Geometry: POINT (0 180)

julia> ArchGDAL.importEPSGA(4326) do sourcecrs_ref
           ArchGDAL.importPROJ4("+proj=cea") do targetcrs_ref
               ArchGDAL.createcoordtrans(sourcecrs_ref, targetcrs_ref) do transform
                   ArchGDAL.transform!(p, transform)
               end
           end
       end
Geometry: POINT (20037508.3427892 0.0)

For my use cases it breaks the generality or reproject when the axis order switches
for EPSG, which worked when I first wrote it! I'm not sure how to handle that now. Probably
all upstream uses should convert EPSG to PROJ before using reproject.

[visr]

I see, that docstring must be updated.

From what I understand, PROJ now follows the ordering as defined by the authority. For EPSG that is typically lat/lon, see https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent.

PROJ also offers proj_normalize_for_visualization, see also https://proj.org/development/quickstart.html. Would you want to such a function? Not sure how/if this is exposed through the GDAL API.

[rafaqz]

proj_normalize_for_visualization could help. If it's accessible through GDAL we could dispatch to using it by defining a geometry type that guaranteed an axis order. That would be useful and let us do generic things without worrying about the differences.

If it's not accessible we could define axis order traits for the formats in GeoFormatTypes. Then we can at least automate swapping them back and forth. Although maybe the CRS string actually needs to be parsed to know the return order for some formats?

However it happens, I think we need a way to buffer users from this when that is useful.

Relevent discussion.
r-spatial/sf#1033

[rafaqz]

Ok so WKT does this as well - but unlike EPSG it needs to be parsed to check the authority to know if it will happen. If you convert EPSG(4326) to WellKnownText a coord is reversed by transform!, but if you convert it to Proj it isn't because the authority isn't included. This is kind of horrific the more I think about it.

[rafaqz]

We can also set OSRAxisMappingStrategy to OAMS_TRADITIONAL_GIS_ORDER. A keyword argument in ArchGDAL methods could give users control over this.

[visr]

We can also set OSRAxisMappingStrategy to OAMS_TRADITIONAL_GIS_ORDER. A keyword argument in ArchGDAL methods could give users control over this.

I see, this is explained in https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order. For the C API, in the ogr_srs_api there are also functions and enums available, search AxisMapping. Those have been wrapped in GDAL.jl, but are not in ArchGDAL.jl.

From your sf link I also found this on the R-spatial blog: https://www.r-spatial.org/r/2020/03/17/wkt.html:

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):

It's a bit of a tough call to make, but right now I'm thinking that in a GDAL Julia package we should follow GDAL default behavior. Like sf, I can understand that other packages don't want to deal with that and force GIS order. So in any case we should make that easy to do.

Ok so WKT does this as well - but unlike EPSG it needs to be parsed to check the authority to know if it will happen. If you convert EPSG(4326) to WellKnownText a coord is reversed by transform!, but if you convert it to Proj it isn't because the authority isn't included.

Yeah I think this is one of the reasons PROJ discourages using proj strings for CRS. Note though that knowing something is an EPSG code is not enough. The rule is just that EPSG decides per projection, hence (from the same R-spatial blog post):

PROJ respects authorities (such as EPSG) for determining whether coordinate pairs refer to longitude-latitude (such as 3857), or latitude-longitude (such as 4326)

[rafaqz]

a GDAL Julia package we should follow GDAL default behavior.

I agree. Some wrapping packages can change that like sf does. ArchGDAL can have a kwarg for it, the question is if it defaults to traditional or mixed order.

When I get around to handing point data in GeoData.jl it will probably use the custom order option - the order will be arbitrary but fixed for any one object depending on the order of the dimensions.