Lunar elevation inconsistent between osgEarth and GDAL

[ravidavi]

Hi all, after a long hiatus I'm getting back into using osgEarth to model the Moon. However, I'm finding an elevation inconsistency between osgEarth and GDAL, specifically gdallocationinfo. For reference, I'm using osgEarth 2.10.2 and GDAL 2.1.3; yes I will upgrade to osgEarth 3.x, but need to confirm historical capability first.

I'm using a Lunar DEM from https://astrogeology.usgs.gov/search/map/moon_lro_lola_dem_118m (download link at top). Warning it's a large 8GB file.

I tile the DEM using: (warning this takes a long time)

osgearth_conv.exe --in driver gdal --in url "Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif" --out driver tms --out format tif --out url "moontiles\tms.xml" --threads 8 --elevation

For testing you can limit the tiling area around (lat,lon) = (0,-90) with --extents -10 -100 10 -80

My corresponding moon.earth file is below. Use any moon image if you don't want to look at a white sphere.

<map name="Moon" type="geocentric" version="2">  
  <options>
    <profile srs="+proj=longlat +R=1737400 +units=m +no_defs"/>
  </options>
  
  <elevation name="terrain" driver="tms">
    <url>moontiles\tms.xml</url>
  </elevation> 
</map>

If I open this with osgearth_viewer.exe --coords moon.earth, then I can manually zoom into (lat,lon) (0,-90). If I do this enough (to ensure loading the most detailed terrain tile), then I observe an elevation of 358.56m (from osgearth_viewer).

In contrast, if I directly check the elevation of the DEM via gdallocationinfo.exe "Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif" -wgs84 -90 0, then I see an elevation of 336m.

So my question is, why the elevation difference? Is this an osgEarth 2.10.2 thing? Am I tiling the DEM wrong, or referencing it incorrectly from my earth file?

Thanks in advance!
Ravi

[ravidavi]

FYI, the output of gdalinfo:

Driver: GTiff/GeoTIFF
Files: Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif
Size is 92160, 46080
Coordinate System is:
PROJCS["SimpleCylindrical Moon",
    GEOGCS["GCS_Moon",
        DATUM["D_Moon",
            SPHEROID["Moon",1737400,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-5458203.076608000323176,2729101.538304000161588)
Pixel Size = (118.450587600000006,-118.450587600000006)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-5458203.077, 2729101.538) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Left  (-5458203.077,-2729101.538) (180d 0' 0.00"E, 90d 0' 0.00"S)
Upper Right ( 5458203.077, 2729101.538) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Right ( 5458203.077,-2729101.538) (180d 0' 0.00"W, 90d 0' 0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=92160x1 Type=Int16, ColorInterp=Gray
  NoData Value=-32768
  Offset: 0,   Scale:0.5

[gwaldron]

Hi Ravi,

I suspect the discrepency you see is down to two things:

1. Either osgEarth 2.x or GDAL 2.1.3 was not properly accounting for the GeoTIFF scale factor (a bug)
2. You need to tell gdallocationinfo to use bilinear interpolation.

I tested this using GDAL 3.10.
QGIS and Global Mapper report the correct value at Lon Lat (-90,0) = ~179m.

Run gdallocationinfo with bilinear interpolation:

set PROJ_IGNORE_CELESTIAL_BODY=1
gdallocationinfo Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.tif -wgs84 -r bilinear -90 0
Report:
  Location: (23040.0000011021P,23040L)
  Band 1:
    Value: 358.499995591526
    Descaled Value: 179.249997795763

The descaled value is what you want. gdalinfo says:

Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
  Min=-18121.000 Max=21428.000
  NoData Value=-32768
  Overviews: 46080x23040, 23040x11520, 11520x5760, 5760x2880, 2880x1440, 1440x720
  Offset: 0,   Scale:0.5

That Scale:0.5 is why you see 358 instead of 179.

Running all this is osgEarth 3.7 with GDAL 3.10, the values line up very closely. A small (probably submeter) error is not unexpected due to the resampling of the heightfields for rendering. Hope this helps!

[gwaldron]

Visualizing the GeoTIFF directly, the values look correct.

Are you getting the 358 with TMS? If so, did you recreate the TMS repo with osgEarth 3.7 or are you still using the one created under osgEarth 2?

[ravidavi]

Aah yeah I was using the one created with osgEarth 2.10.2. I'm trying to recreate it with 3.7 but osgearth_conv is now taking 10x longer than it did with 2.10.2. Not sure why (maybe a vcpkg thing?), but I'll start another discussion about that.

BTW how are you able to visualize the geotiff directly? Isn't an 8GB geotiff too large to use in osgEarth without tiling first?

[gwaldron]

It's best to preprocess the GeoTIFF before using it, by reprojecting, adding internal tiles, and adding overviews.
Examples below (from memory so forgive any typos)...

Reproject it

gdalwarp -t_srs "+proj=longlat +R=1737400 +units=m +no_defs" input.tif temp1.tif

Add internal tiling:
(edit: you can add compression, but it slows down the reading a lot)

gdal_translate -of GTiff -co TILED=YES -co BIGTIFF=YES -co NUM_THREADS=16 temp1.tif output.tif

Add overviews:

gdaladdo -r average --config GDAL_NUM_THREADS 16 --config COMPRESS_OVERVIEW DEFLATE output.tif 2 4 8 16 32 64

It is still slow but quite a bit faster. You can also enable caching in osgEarth..

set OSGEARTH_CACHE_PATH=cache

If you would rather tile it you can, and I would recommend using MBTilesElevation as your output format instead of TMSElevation. It's just less cumbersome having a single file.

There are more details at https://docs.osgearth.org/en/latest/data.html#optimizing-a-geotiff

[ravidavi]

All these years I thought I had to tile a very large geotiff, because loading the 8GB raw file would not be possible! I guess I was assuming the whole file needed to be loaded; lesson learned thank you!

Out of curiosity, in the osgearth_imgui terrain panel, what is the Z in the first output line Lat ... Lon ... Z ...? I notice that it matches the HAE elevation when zoomed in at most of the globe, but it differs from HAE near the poles (no matter how much I zoom in). Is that difference expected?

[gwaldron]

Well, with the combination of -co TILED=YES and gdaladdo you effectively ARE tiling it, just internally to the GeoTIFF format :)

You will notice though, that even with tiling, it's a bit slow reading the DEM through GDAL. That is because the current algorithm uses point-by-point interpolation to transform the source data to the triangle mesh, and this is quite slow. We are working on a new approach to speed this up.

As for the terrain panel:

As you know the "HAE Elevation" queries the source data DEM directly, and reports the resolution of the resulting value. (At the moment, that resolution will be wrong for the Moon. It's assuming an Earth radius. And the MSL readout only applies to the Earth.)

That "Z" readout queries the OpenGL geometry (the actual triangles in the rendered mesh). The mesh is a lo-res approximation of the DEM at various levels of detail, and thus the accuracy increases the closer you zoom and the more dense the triangle mesh gets.

Hope this helps :)

[XenonofArcticus]

Ravi, if you are working on Lunar stuff, I have done some work in this field and would be happy to chat about it.

…

On Mon, Dec 23, 2024 at 1:20 PM Glenn ***@***.***> wrote: I see 4 pixels that touch exactly at the location (-90,0). Their values are 187.5 188.5 173.0 168.0 Yielding a bilinear interpolation = 179.25. Perhaps QuickMap is applying a half-pixel shift to its sampling as one might commonly do for imagery (Area v. Point), and landing directly on the 168. Or perhaps the other 3 softwares are not applying it when they should. Hard to say without ground truth! — Reply to this email directly, view it on GitHub <#2672 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAPFE3HS44QSSOTJYY32NTD2HBV7NAVCNFSM6AAAAABUCDLNZCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCNRVGMYTQOA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[ravidavi]

Hey Chris, yes I'd love to chat about it after the holidays! I'll email you and we can set up a time.