<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">Dear all,<div><br></div><div>Question: how do I define my source raster instance, to avoid using the advance `-ct` operator in gdalwarp in GDAL3 during reprojection?</div><div><br></div><div>I've been reading <a href="https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn#Axisorderissues">https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn#Axisorderissues</a> and <a href="https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html">https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html</a> and I've tried numerous things and I really wished I could have found the info by myself. </div><div>But I'm a bit lost.</div><div><br></div><div>I've a .h5 file with data originating from the Dutch weather radar. I've an <i>old fashion</i> 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.</div><div><br></div><div>But now, using GDAL3, I cannot get the gdalwarp reprojection to give me anything that comes close to what was possible in GDAL2.</div><div><br></div><div>In GDAL2 the following works:</div><div><pre style="box-sizing:border-box;font-family:SFMono-Regular,Consolas,"Liberation Mono",Menlo,monospace;font-size:11.9px;margin-top:0px;padding:16px;overflow:auto;line-height:1.45;border-top-left-radius:6px;border-top-right-radius:6px;border-bottom-right-radius:6px;border-bottom-left-radius:6px;color:rgb(36,41,46);font-variant-ligatures:normal;margin-bottom:0px"><code style="box-sizing:border-box;font-family:SFMono-Regular,Consolas,"Liberation Mono",Menlo,monospace;font-size:11.9px;padding:0px;margin:0px;border-top-left-radius:6px;border-top-right-radius:6px;border-bottom-right-radius:6px;border-bottom-left-radius:6px;word-break:normal;border:0px;display:inline;overflow:visible;line-height:inherit">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  
</code></pre><br class="gmail-Apple-interchange-newline" style="color:rgb(0,0,0)"></div><div><pre style="box-sizing:border-box;font-family:SFMono-Regular,Consolas,"Liberation Mono",Menlo,monospace;font-size:11.9px;margin-top:0px;padding:16px;overflow:auto;line-height:1.45;border-top-left-radius:6px;border-top-right-radius:6px;border-bottom-right-radius:6px;border-bottom-left-radius:6px;color:rgb(36,41,46);font-variant-ligatures:normal;margin-bottom:0px"><code style="box-sizing:border-box;font-family:SFMono-Regular,Consolas,"Liberation Mono",Menlo,monospace;font-size:11.9px;padding:0px;margin:0px;border-top-left-radius:6px;border-top-right-radius:6px;border-bottom-right-radius:6px;border-bottom-left-radius:6px;word-break:normal;border:0px;display:inline;overflow:visible;line-height:inherit">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</code></pre></div><div><br></div><div>But now with GDAL3 using the same command I get:</div><div><br></div><div><p style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo;color:rgb(0,0,0)"><span style="font-variant-ligatures:no-common-ligatures">ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body</span></p></div><div><br></div><div>And using the `-ct` operator with GDAL3 I get this:</div><div><br></div><div><p style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo;color:rgb(0,0,0)"><span style="font-variant-ligatures:no-common-ligatures">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</span></p>
<p style="margin:0px;font-stretch:normal;font-size:11px;line-height:normal;font-family:Menlo;color:rgb(0,0,0)"><span style="font-variant-ligatures:no-common-ligatures">0 -3650 0</span></p></div><div><br></div><div>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. </div><div><br></div><div>I've prepared all this as well in this GitHub <a href="https://github.com/mattijn/KNMI-radar-PyGMT/blob/main/reproject%20knmi%20radar%20using%20rasterio.ipynb">Notebook</a> 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.</div><div><br></div><div>Any help or advice is much appreciated!</div><div><br></div><div>Kind regards,</div><div>Mattijn van Hoek</div></div></div></div></div></div></div></div></div></div>