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.