ArchGDAL transforming CRS incorrectly (probably wrong long/lat order)

Using GeoDataFrames.jl I’m trying to reproject to Google Pseudo-Mercator (EPSG 3857) projection, but get wrong coordinates.

In the original GeoJSON, coords are long/lat (see screenshot below):

                    geometry
0  POINT (69.19970 41.36596)
1  POINT (69.62239 42.41753)
2  POINT (74.57231 42.90830)
3  POINT (76.99721 43.24954)

Julia code:

using GeoFormatTypes; const GFT=GeoFormatTypes
using GeoDataFrames; const GDF=GeoDataFrames

cities = GDF.read("/tmp/cities.geojson")
reproject(cities.geometry, GFT.EPSG(4326), GFT.EPSG(3857))
GDF.write("/tmp/cities-julia.geojson", cities)
cities.geometry

Output:

4-element Vector{ArchGDAL.IGeometry{ArchGDAL.wkbPoint}}:
 Geometry: POINT (4604838.04865289 10813102.9330808)
 Geometry: POINT (4721897.84030841 10946911.4321356)
 Geometry: POINT (4776530.55208298 12750804.4839833)
 Geometry: POINT (4814516.76984332 13852710.5883733)

Python code:

import geopandas as gpd
cities = gpd.read_file('/tmp/cities.geojson').to_crs(3857)
cities.to_file('/tmp/cities-python.geojson')
cities

Output:

                          geometry
0  POINT (7703275.478 5066472.054)
1  POINT (7750329.114 5223730.064)
2  POINT (8301351.354 5298025.179)
3  POINT (8571290.321 5350031.840)

Plotting it on the map, Python reprojected points are showing, but Julia’s don’t appear at all:

It looks like ArchGDAL thinks coords are in lat/lon order, while they’re actually lon/lat. But the docs on ArchGDAL and GeoFormatTypes are very cryptic. GFT docs say nothing of axis order, and I even don’t know if there’s an ID for EPSG 4326 with long/lat sequence.

In ArchGDAL there’s a method .importEPSG ([can’t add more than 2 links per post] on what? on the whole ArchGDAL module?), but it’s a total mystery what it does. At least, calling ArchGDAL.importEPSG(4326) (default sequence is long/lat, according to the docs) changed nothing, the output coords are the same as in the first example.

This won’t work either:

ArchGDAL.reproject(cities2.geometry, ArchGDAL.importEPSG(4326), ArchGDAL.importEPSG(3857))

What do I do?

This is such a simple task, yet googling produces only one answer with the simplest example [can’t add more than 2 links per post], copied from the GeoDataFrames docs.

Ok. So GFT just provides wrappers for CRS so it can be shared between packages without re-parsing. It doesn’t really do anything with it.

GeoDataFrames.jl uses gdal via ArchGDAL to both read and reproject data.

It seems to me your data is being read as x/y rather than y/x at some point in the chain, as its the default. Most likely this is a bug in ArchGDAL.jl.

1 Like

Have you tried the order keyword argument to ArchGDAL.importEPSG? As far as I understand the docstring setting the order to trad could swap the axis to lon/lat.

2 Likes

Indeed this looks like an axis order issue. The reproject function that you are using has a keyword argument order=:trad that you can use. The issue is that EPSG defines the argument order to be lat/lon, which is not the case here. And the axis order is taken based on the authority (in this case EPSG) by default. See also some more details here: FAQ — PROJ 9.1.0 documentation.

3 Likes

Tried this, it didn’t change anything.

But I’m not sure how to use it – should I just call it, or pass the return value somewhere? Passing it to .reproject causes MethodError.

Tried this, but the result is still the same:

:compliant is the default, can you pass order=:trad?

1 Like

It did the trick! Thanks! Hmmm… happy it worked, but sad I have to do such tricks.

I’m still perplexed how cryptic the docs are. No info on return types. Also this example crashed with MethodError, couldn’t find how to call it on a vector.


cities2 = GDF.read("/tmp/cities.geojson")
sc = AG.importEPSG(4326)
tc = AG.importEPSG(3857)

AG.createcoordtrans(sc, tc) do transform
    AG.transform!.(cities2.geometry, transform)
end

It’s best to make github issues for any problems you find.

The docs can always be improved, but there are only a few people using small amounts of our spare time working on these packages. We all have other jobs.

If you feel strongly about this you may like to get involved and for a start improve the docs.

As for the issue, I feel there is actually a bug in ArchGDAL as we treat the order as if it’s always x/y/z when in your case its y/x. Fixing that will fix plotting and other things.

I would prefer to have :trad as the default (which would fix your problem). CRS dependent coordinate order is an implementation nightmare and IMHO should never have been made the standard (in gdal, this is not a julia thing). But others feel otherwise I’m sure.

2 Likes

I fully agree that the order should always be lon,lat and not lat,lon but this is a PROJ thing that can be avoided in GDAL if OAMS_TRADITIONAL_GIS_ORDER is used in OSRSetAxisMappingStrategy()

We have this in GMT (C lib) and AFAIK that is also the default in GMT.jl when it uses PROJ directly.

#if GDAL_VERSION_MAJOR >= 3
	OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER);		/* Set the data axis to CRS axis mapping strategy. */
#endif
1 Like