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

**URL:** <https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098>\
**Category:** New to Julia\
**Tags:** geo, geodesy\
**Created:** [November 11, 2022, 6:41am UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098 "2022-11-11T06:41:27Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![culebron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/culebron/32/44276_2.png) [@culebron](https://discourse.julialang.org/u/culebron)\
**Post date:** [November 11, 2022, 6:41am UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/1 "2022-11-11T06:41:27Z")

</div>

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

```julia
                    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:

```julia
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:

```julia
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:

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

```

Output:

```julia
                          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:

 ![изображение](https://global.discourse-cdn.com/julialang/original/3X/d/5/d50c10e6776541ff394d601014e2e7c37244f12b.jpeg)

It looks like ArchGDAL thinks coords are in lat/lon order, while they’re actually lon/lat. But [the docs on ArchGDAL](https://yeesian.com/ArchGDAL.jl/dev/reference/#ArchGDAL.importEPSG-Tuple%7BInteger%7D) and GeoFormatTypes are very cryptic. [GFT docs](https://juliageo.org/GeoFormatTypes.jl/stable/#GeoFormatTypes.EPSG) 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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [November 11, 2022, 8:22am UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/2 "2022-11-11T08:22:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [November 11, 2022, 8:33am UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/3 "2022-11-11T08:33:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [November 11, 2022, 10:24am UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/4 "2022-11-11T10:24:41Z")

</div>

Indeed this looks like an axis order issue. The [reproject function](https://yeesian.com/ArchGDAL.jl/stable/reference/#ArchGDAL.reproject-Union%7BTuple%7BT%7D,%20Tuple%7BT,%20GeoFormatTypes.GeoFormat,%20Nothing%7D%7D%20where%20T) 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](https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent).

---

<div class="post-metadata">

**Author:** ![culebron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/culebron/32/44276_2.png) [@culebron](https://discourse.julialang.org/u/culebron)\
**Post date:** [November 11, 2022, 12:09pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/5 "2022-11-11T12:09:50Z")

</div>

> [@Fliks](#):
>
> 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.

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`.

> [@visr](#):
>
> Indeed this looks like an axis order issue. The [reproject function](https://yeesian.com/ArchGDAL.jl/stable/reference/#ArchGDAL.reproject-Union%7BTuple%7BT%7D,%20Tuple%7BT,%20GeoFormatTypes.GeoFormat,%20Nothing%7D%7D%20where%20T) that you are using has a keyword argument `order=:trad` that you can use.

Tried this, but the result is still the same:

 ![изображение](https://global.discourse-cdn.com/julialang/original/3X/7/6/76ea7487b9786c5c5d044c14eb9bb2b20a4b1c61.png)

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [November 11, 2022, 12:19pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/6 "2022-11-11T12:19:08Z")

</div>

> [@culebron](#):
>
> Tried this, but the result is still the same

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

---

<div class="post-metadata">

**Author:** ![culebron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/culebron/32/44276_2.png) [@culebron](https://discourse.julialang.org/u/culebron)\
**Post date:** [November 11, 2022, 12:20pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/7 "2022-11-11T12:20:40Z")

</div>

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

 ![изображение](https://global.discourse-cdn.com/julialang/original/3X/a/1/a1b4a1dcc9a869fa176e3c0fb5bb425ea23650e7.png)

---

<div class="post-metadata">

**Author:** ![culebron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/culebron/32/44276_2.png) [@culebron](https://discourse.julialang.org/u/culebron)\
**Post date:** [November 11, 2022, 12:26pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/8 "2022-11-11T12:26:45Z")

</div>

I’m still perplexed how cryptic the docs are. No info on return types. Also [this example](https://discourse.julialang.org/t/archgdal-jl-or-gdal-jl-not-sure-undesired-behaviour-when-reprojecting-coordinates-to-epsg-4326/59415) crashed with MethodError, couldn’t find how to call it on a vector.

```julia

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

```

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [November 11, 2022, 5:16pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/9 "2022-11-11T17:16:43Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [November 11, 2022, 5:36pm UTC](https://discourse.julialang.org/t/archgdal-transforming-crs-incorrectly-probably-wrong-long-lat-order/90098/10 "2022-11-11T17:36:59Z")

</div>

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.

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

```