[gdal-dev] OGRCreateCoordinateTransformation problem on GDAL 3.0.3 and PROJ 6.3.0

[Edzer Pebesma]

On 1/23/20 5:58 PM, Even Rouault wrote:
>> Thanks! Yes, with SRS->importFromUserInput("EPSG:3857"); etc this works
>> fine; I'm OK with the warning, I just don't see how the subsequent error
>> relates to syntax that is deprecated.
> 
> I'm looking at this, but haven't yet an explanation. This is very subtle,
> and seems to be linked specifically to EPSG:3857

I don't think so: I see the same problem if I use +init=epsg:28992
rather than +init=epsg:3875.

Your advice of improving this C++ program is well meant, but fixing this
C++ program is not my problem, the problem is how we use GDAL's SRS
interface in R packages. I saw the error occuring at pretty random
occasions, but always errors _after_ one or more +init=epsg:xxxx strings
had been used to initialize an SRS.

> 
> I've modified your reproducer as the following:
> 
> int main() {
> 
>    OGRSpatialReference *aSRS = new OGRSpatialReference;
> 
>    aSRS->importFromWkt("GEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]");
> 
>    OGRSpatialReference *bSRS = new OGRSpatialReference;
> #ifdef WORKING
>    OGRSpatialReference srs;
>    srs.importFromProj4("+init=epsg:3857");
>    char* wktb = NULL;
>    srs.exportToPrettyWkt(&wktb);
> #else
>    bSRS->importFromProj4("+init=epsg:3857");
>    const char* wktb = "PROJCRS[\"WGS 84 / Pseudo-Mercator\",BASEGEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"unnamed\",METHOD[\"Popular Visualisation Pseudo Mercator\",ID[\"EPSG\",1024]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]";
> #endif
>    bSRS->importFromWkt((const char *) wktb);
> 
>    OGRCoordinateTransformation *ct =
> OGRCreateCoordinateTransformation(aSRS, bSRS);
>    if (ct == NULL) {
>         printf("ct NULL\n");
>         exit(1);
>    }
>    exit(0);
> }
> 
> So the error is linked to having importFromProj4() and then importFromWkt()
> This doesn't make sense as the later should cancel the former, so there's some
> side effect of importFromProj4() that has later consequences.
> If you define WORKING, which will do the importFromProj4() + exportToPrettyWkt() 
> in a temporary object, and import the resulting WKT in the final bSRS, then
> it works. Doesn't make sense either but could be a workaround
> 
> But Sean's advice is definitely the way forward . Use importFromEPSG(XXXX) or
> SetFromUserInput("EPSG:XXXX")
> And SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) if you don't
> want EPSG compliant axis definitions.
> 
> Even
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 3110 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200123/1f133c4f/attachment.key>