[gdal-dev] Troubles in coord conversions with GDAL+proj6

[Joaquim Manuel Freire Luís]

Not from the same command as before because of that Matlab-filtering-stderr issue that I mentioned but this should be equivalent. It's from Mirone (compiled as true stand alone exe) and that uses calls to the same ogrproj code. Not exactly what you expected, I think. 

-- read a geotiff file

GDAL: GDALOpen(C:\v\lixo.tiff, this=0189FE78) succeeds as GTiff.
ERROR 1: PROJ: pj_open_lib(proj.db): call fopen(C:\SVN\mirone30\data\proj_lib\proj.db) - succeeded
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDAL_CACHEMAX = 102 MB
GDAL: GDALClose(C:\v\lixo.tiff, this=0189FE78)
GDAL: GDALOpen(C:\v\lixo.tiff, this=0189FE78) succeeds as GTiff.
GDAL: GDALClose(C:\v\lixo.tiff, this=0189FE78)

The ogrproj call should start here

ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=WGS84
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=WGS84
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
ERROR 1: PROJ: xy_in unit: Radian
ERROR 1: PROJ: xy_out unit: Degree
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=WGS84
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=WGS84
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
ERROR 1: PROJ: xy_in unit: Radian
ERROR 1: PROJ: xy_out unit: Degree
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=WGS84
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
ERROR 1: PROJ: pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
ERROR 1: PROJ: pj_ellipsoid - final:    ellps=GRS80
...


|>-----Original Message-----
|>From: Even Rouault <even.rouault at spatialys.com>
|>Sent: Sunday, February 24, 2019 1:29 PM
|>To: gdal-dev at lists.osgeo.org
|>Cc: Joaquim Manuel Freire Luís <jluis at ualg.pt>
|>Subject: Re: [gdal-dev] Troubles in coord conversions with GDAL+proj6
|>
|>Joaquim,
|>
|>>
|>> I'm having troubles migrating some of my Matlab mex to recent GDAL
|>> (built yesterday from git). With older GDAL versions all went well.
|>> For example, this is correct
|>> >> prj.SrcProjSRS = '+proj=utm +zone=29 +datum=WGS84 +units=m
|>> >> +no_defs'; xy = ogrproj([0 0], prj)
|>>
|>> xy =
|>>
|>>  -13.488743884387196                   0
|>>
|>> But with newest version I'm having two problems
|>>
|>> >> prj.SrcProjSRS = '+proj=utm +zone=29 +datum=WGS84 +units=m
|>> >> +no_defs'; xy = ogrproj([0 0], prj)
|>>
|>> xy =
|>>
|>>   -0.000121660050290 -13.488743885487231
|>>
|>> One is that coordinates are not correct (lat should be 0.0)
|>
|>That one is really odd. It looks like some datum shift would be applied, which
|>doesn't make sense here since your target SRS uses WGS84 too. Can you run
|>your code with PROJ_DEBUG=5 and CPL_DEBUG=ON environment variables set
|>?
|>
|>For correct results, you should see something like:
|>OGRCT: proj=pipeline step inv proj=utm zone=29 ellps=WGS84 step
|>proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap order=2,1 (Inverse
|>of UTM zone 29N)
|>
|>
|>> and the other is
|>> that in this later case it outputs lat,long instead of lon,lat
|>
|>That one is expected. SetWellKnownGeogCS( "WGS84" ) now returns a SRS
|>with lat, lon order. You can force "GIS" friendly order by calling
|>oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
|>
|>Even
|>
|>--
|>Spatialys - Geospatial professional services http://www.spatialys.com