[gdal-dev] Transformation works in GDAL2 but cause problems in GDAL3

[mattijn]

Dear all,

Question: how do I define my source raster instance, to avoid using the
advance `-ct` operator in gdalwarp in GDAL3 during reprojection?

I've been reading
https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn#Axisorderissues
 and
https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html
and I've tried numerous things and I really wished I could have found the
info by myself.
But I'm a bit lost.

I've a .h5 file with data originating from the Dutch weather radar.
I've an *old
fashion* PROJ.4 definition and an Affine instance. Based on this
information I was, using GDAL2, able to create a registered-raster instance
or geotiff. And then using gdalwarp I could reproject it to a more
common epsg:28992 or epsg:4326 projection.

But now, using GDAL3, I cannot get the gdalwarp reprojection to give me
anything that comes close to what was possible in GDAL2.

In GDAL2 the following works:

echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0
+y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs
EPSG:4326 --debug on


OGRCT: Source: +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0
+y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0
+y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs
0 55.9735620711571 0


But now with GDAL3 using the same command I get:

ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do not
belong to the same celestial body

And using the `-ct` operator with GDAL3 I get this:

echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0
+lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 -ct
'+proj=axisswap +order=2,1 +step' --debug on

0 -3650 0

All in all, I like to know how I could create or define my source raster
instance operable for GDAL3 in such a way, it is able to reproject it
correctly without any advanced reprojection operators.

I've prepared all this as well in this GitHub Notebook
<https://github.com/mattijn/KNMI-radar-PyGMT/blob/main/reproject%20knmi%20radar%20using%20rasterio.ipynb>
in Python, that works with GDAL2 installed, but fails similar as above with
GDAL3. In this notebook is referred to a file that is included in the repo
of the notebook.

Any help or advice is much appreciated!

Kind regards,
Mattijn van Hoek
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201202/9a2a0871/attachment-0001.html>