Dem size is not computed correctly for extraterrestrial celestial bodies

[scpeters]

While debugging #2828 (a hard failure of gazebo when loading certain Dem files), I noticed the following error message from the reported UNIT_Dem_TEST failure:

[ RUN      ] DemTest.NegDem
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["Moon2000_npole",BASEGEOGCRS["GCS_Moon",DATUM["Moon_2000",ELLIPSOID["Moon_2000_IAU_IAG",1737400,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",90,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["easting",south,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1]]]'
 to `EPSG:4326'

Based on the first error message "Source and target ellipsoid do not belong to the same celestial body" and noting that EPSG:4326 refers to the WGS84 model of earth and that the PROJCRS string for dem_neg.tif used in this test appears to refer to the moon, there are clearly issues with the coordinate transformations used in our Dem class.

The following is the output of gdalinfo dem_neg.tif -mm:

Driver: GTiff/GeoTIFF
Files: ../test/media/materials/textures/dem_neg.tif
Size is 33, 33
Coordinate System is:
PROJCS["Moon2000_npole",
    GEOGCS["GCS_Moon",
        DATUM["Moon_2000",
            SPHEROID["Moon_2000_IAU_IAG",1737400,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-53750.000000000000000,-51050.000000000000000)
Pixel Size = (2.426609848484849,-2.426609848484849)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -53750.000,  -51050.000) ( 46d28'32.90"W, 87d33'20.67"N)
Lower Left  (  -53750.000,  -51130.078) ( 46d25'51.46"W, 87d33'14.12"N)
Upper Right (  -53669.922,  -51050.000) ( 46d25'59.33"W, 87d33'27.55"N)
Lower Right (  -53669.922,  -51130.078) ( 46d23'17.88"W, 87d33'21.00"N)
Center      (  -53709.961,  -51090.039) ( 46d25'55.39"W, 87d33'20.84"N)
Band 1 Block=33x33 Type=Float32, ColorInterp=Gray
    Computed Min/Max=-212.296,-205.440

Based on its 33x33 image size and Pixel Size = (2.426609848484849,-2.426609848484849), I would expect the world height and width of this Dem to be 33 * 2.426609848484849 ~= 80.078125, but Dem_TEST expects 293.51068 for the height and width.

[scpeters]

A workaround is to use the <size> element to manually specify the heightmap length and width

[scpeters]

The Dem::Load function currently calls GetGeoReference for different corners of the Dem and computes a distance between those points using the SphericalCoordinates class. The GetGeoReference function tries to convert everything into WGS84 to match with the assumptions of the SphericalCoordinates class, which can work for earth, but I believe this operation is failing for Dem's based on other celestial bodies.

[scpeters]

the test/worlds/dem_neg.world file currently seg-faults without a <size> element, so I've added one in scpeters@ba64bc7

[scpeters]

the test/worlds/dem_neg.world file currently seg-faults without a <size> element, so I've added one in scpeters@ba64bc7

current console output with this change:

$ gzserver --verbose test/worlds/dem_neg.world 
Gazebo multi-robot simulator, version 11.10.2
Copyright (C) 2012 Open Source Robotics Foundation.
Released under the Apache 2 License.
http://gazebosim.org

[Msg] Waiting for master.
[Msg] Connected to gazebo master @ http://127.0.0.1:11345
[Msg] Publicized address: 172.27.1.28
[Msg] Loading world file [/home/scpeters/clone/gazebo/test/worlds/dem_neg.world]
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["Moon2000_npole",BASEGEOGCRS["GCS_Moon",DATUM["Moon_2000",ELLIPSOID["Moon_2000_IAU_IAG",1737400,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",90,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["easting",south,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1]]]' to `EPSG:4326'
[Err] [Dem.cc:236] EXCEPTION: Unable to transform terrain coordinate system to WGS84 for coordinates (0,0)

[Wrn] [Dem.cc:133] Failed to automatically compute DEM size. Please use the <size> element to manually set DEM size.
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["Moon2000_npole",BASEGEOGCRS["GCS_Moon",DATUM["Moon_2000",ELLIPSOID["Moon_2000_IAU_IAG",1737400,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference_Meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",90,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["easting",south,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1]]]' to `EPSG:4326'
[Err] [Dem.cc:236] EXCEPTION: Unable to transform terrain coordinate system to WGS84 for coordinates (0,0)

[Wrn] [HeightmapShape.cc:125] DEM coordinate transformation error. SphericalCoordiantes and GpsSensor may not function properly.

[adityapande-1995]

The error is thrown by this particular call in Dem.cc : https://github.com/osrf/gazebo/blob/72d47eb0f08e24d8a5c570d26c6e360d97ab9242/gazebo/common/Dem.cc#L233-L238

which is implemented here in GDAL : https://github.com/OSGeo/gdal/blob/4a5f171a2beb6c29fd5dd8b4044cd2539fce6b22/ogr/ogrct.cpp#L1514-L1529

[adityapande-1995]

This problem has been addressed here : gazebosim/gz-sim#1556
We can add new surfaces : MOON_SCS and CUSTOM_SURFACE to indicate celestial bodies.

[adityapande-1995]

This issue has been fixed in Gazebo garden and needs to be backported to gazebo classic.