[gdal-dev] warp issue: Source and target ellipsoid do not belong to the same celestial body

[Even Rouault]

Le 22/01/2024 à 17:33, Wilco K via gdal-dev a écrit :
> Hi,
>
> these 2 lines work with GDAL 2.2.3, but not with GDAL 3.5.0 anymore.
>
> gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 
> +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -a_ullr 0.0, 
> -3649.9792, 700.000, -4414.9792 
> "HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
> "test-translate.tif"
>
> gdalwarp -t_srs EPSG:4326 -dstnodata 65535 -of GTiff 
> "test-translate.tif" "test-warp.tif"
>
>
> Error:
> ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do 
> not belong to the same celestial body
> ERROR 6: Cannot find coordinate operations from 
> `PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6378.14,298.183263207102,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference 
> meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["Polar 
> Stereographic (variant B)",METHOD["Polar Stereographic (variant 
> B)",ID["EPSG",9829]],PARAMETER["Latitude of standard 
> parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude 
> of 
> origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],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)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1]]]' 
> to `EPSG:4326'
>
>
> But when the +a and +b in the gdal_translate call are changed, 
> gdalwarp does work:
> gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 
> +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0" -a_ullr 0 
> -3649999.11191775 700000.903671186 -4415003.88199764 
> "HDF5:\"RAD_NL25_RAC_03H_202401220800.h5\"://image1/image_data" 
> "test-translate.tif"
>
> What is the problem?

The value of the +a and +b parameters must be in metres, not in km. So 
you add to cheat also on the geotransform. But if datum transformations 
had to be done the fact that you transform between an ellipsoid of ~ 
6000 meters to one of ~ 6000 km wouldn't work well... It probably sort 
of worked because datum transformation was skipped, but the new 
behaviour which checks that the shape of the source and target ellipsoid 
isn't too different is definitely saner and will avoid potential 
reprojection errors

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240122/e276e0e0/attachment.htm>