Allow use of elevation data (rather than radius) in a shapemodel

[rvwagner]

Description
Add the ability to use DEM cubes with terrain expressed as elevation relative to the planetary radius as shapemodels, instead of requiring terrain to be expressed as radius. Requiring a radius means that 32-bit files on the Moon will be quantized down to steps of ~10-20cm, and likely worse for larger planets like Mars, which is significantly worse than the level of detail in LROC NAC or HIRiSE DEMs (see #4603).

An alternate solution might be to allow the use of the Base/Multiplier terms in 32-bit files, but that feels like it would have even further-reaching repercussions.

[jessemapel]

I have some questions about how to properly handle this:

1. What datum do we take the elevation values above?
2. How do we specify the datum in the DEM file?

It sounds like you're recommending a spherical datum. We could use something like:

DatumRadius = XXXXX <m>

in the DEM label to specify that. For more complex datums I'm not sure what we'd use. I don't have a good grasp on what datums are commonly used for elevation DEMs. @thareUSGS @dpmayerUSGS any inputs here?

Allowing base and multiplier in 32 bit files sounds reasonable to me as the cube I/O automatically convert everything to a double for calculations.

[rbeyer]

in the DEM label to specify that. For more complex datums I'm not sure what we'd use. I don't have a good grasp on what datums are commonly used for elevation DEMs. @thareUSGS @dpmayerUSGS any inputs here?

Commonly? Ha ha, that's a good one for a Monday morning. This is planetary science: there are no "common" datums (unless you intentionally limit yourself to one planet, depending on the planet--but then that's not a general solution), and the unfortunate short answer is if we're going to support this, then we need to be able to support a mechanism for a user to describe a completely arbitrary datum. No short-cuts for us.

I'm assuming this is for spiceinit right? Does this come up elsewhere in the codebase for a user to specify that? Understanding the impact might also be a lever arm as we think about how to address this.

The "beauty" of the give-me-radius-values approach is that it side-steps the whole how-do-you-describe-a-datum problem. However, the down-side as @rvwagner points out, is that we have a precision problem. Attempting to implement some kind of Base/Multiplier scheme would allow us to keep just taking radii, and would allow more precision, but has the down-side of being non-standard and adding complexity to the input specification (and potentially adding additional confusing parameters to spiceinit or required labels in any DEM provided). Maybe there's some other solution?

[jessemapel]

There are some applications that take a surface model as input, slpmap is one example, and they use a different interface than the camera models and spiceinit use.

The camera models use the ShapeModel class, specifically DemShape for DEMs, to handle this stuff.

Things like slpmap just open up the cube and it's assumed that the values are radii.

[dpmayerUSGS]

I agree with everything @rbeyer said and I too would like to know more about the use case.

How does ASP's mapproject handle this issue? A user has created a high-precision DTM in ASP from a pair of LROC images and now they want to orthorectify the images onto the DTM. In ISIS, camp2map defaults to projecting spiceinit'd images onto a global shape model for the target body (whose pixel values = radius from center of mass to surface).

[rbeyer]

@dpmayerUSGS that is a good question. I do not have the answer to that on the tip of my tongue, but @oleg-alexandrov might.

[KrisBecker]

The Mercury DEM ($ISISDATA/base/dems/MSGR_DEM_USG_EQ_I_V02_prep.cub) uses elevation values with a spherical reference radius (2439400.0 m) specified in the Base keyword and data is rendered as an Equirectangular map projected cube. The pixel type is two-byte signed integers which bins elevation in ~65532 unique elevation values. For this case where the the DEM has a pixel resolution of 665.2 <meters/pixel>, this scenario works just fine. Here is a portion of the labels the describes the DEM in terms of stored elevation:

Object = IsisCube
  Object = Core
    StartByte   = 65537
    Format      = Tile
    TileSamples = 562
    TileLines   = 823

    Group = Dimensions
      Samples = 23042
      Lines   = 11522
      Bands   = 1
    End_Group

    Group = Pixels
      Type       = SignedWord
      ByteOrder  = Lsb
      Base       = 2439400.0
      Multiplier = 1.0
    End_Group
  End_Object

  Group = BandBin
    Name   = MedianRadius
    Number = 5
    Center = 5
    Width  = 5
  End_Group

  Group = Mapping
    ProjectionName       = Equirectangular
    CenterLongitude      = 180.0
    TargetName           = Mercury
    EquatorialRadius     = 2439400.0
    PolarRadius          = 2439400.0
    LatitudeType         = Planetocentric
    LongitudeDirection   = PositiveEast
    LongitudeDomain      = 360
    MinimumLatitude      = -90.0
    MaximumLatitude      = 90.0
    MinimumLongitude     = 0.0
    MaximumLongitude     = 360.00
    UpperLeftCornerX     = -7664266.3623196 <meters>
    UpperLeftCornerY     = 3832465.8027362 <meters>
    PixelResolution      = 665.24315270546 <meters/pixel>
    Scale                = 64.0 <pixels/degree>
    CenterLatitude       = 0.0
    CenterLatitudeRadius = 2439400.0
  End_Group
End_Object

Note the Base keyword specifies the reference radius with Multiplier = 1 - the actual radius is simply (Base + DN).

So such DEMs do exist and can provide the means to store elevation, but this will not directly resolve the OP situation. There are at least two ways to expand this technique to address large bodies: 1) honor the Core Base and Multiplier values in floating point cubes, and 2) add double precision as a supported cube pixel type.

Note that option 1 should be rather easy to implement with little or no complications. Option 2 creates a file size/storage problem that impacts a lot of things, but I still think it would be a good enhancement to ISIS. If size becomes too much of a problem, we can start discussions on how to break up DEMs into smaller files and figure out the best way to support that in ISIS (e.g., a list of DEMs with prioritized ray tracing ).

Also wanted to mention that we routinely swap out DTMs without using spiceinit in our processing methods for various reasons.

[oleg-alexandrov]

In ASP, and likely in every geospatial tool, a DEM is stored as a a GeoTiff with each elevation being a 32-bit floating point number relative to a datum specified in the geoheader. ASP stores point clouds either as 32 bit float coordinates but relative to a (somewhat arbitrarily chosen) "representative nearby point", or using 64 bits relative to the planet center.

When mappoject is invoked, the 32-bit height above datum is converted to a double-precision ECEF value using the datum in the geoheader.

The solution for ISIS is likely to choose a datum and record what datum was chosen.

[dpmayerUSGS]

@oleg-alexandrov Thanks for explaining how ASP does things!

An analogous way to handle things in ISIS could be to use the A/B/C radius values in the mapping group as the parameters for the datum, but this would require an interface change to get programs like spiceinit to deviate from the default behavior of assuming pixels = radius values. There are also situations where the radii used to describe the horizontal projection may be different from the radii defining the datum that pixels are relative to. (This is a general problem; I've been bit by this on 1 or 2 occasions when not being careful handling HRSC DTMs in ASP).

I think I like @KrisBecker 's approach of respecting the Base and Multiplier keywords.

[jessemapel]

@oleg-alexandrov 's comment makes me wonder if this has been handled by geotiffs and using something like GDAL to support geotiff DEMs could resolve this. We could support radius or sphere+elevation Cubes and then use GDAL to support other DEM types.

[rvwagner]

Adding Base/Multiplier support for 32-bit data would certainly solve the problem for my use case (lunar data, LROC NAC DEMs), though it would leave a bit of a gotcha for untrained users, as the most obvious way to turn a cube from elevation to radius (fx or algebra) wouldn't produce the best results, and the best way (editlab) is generally not a function most users should ever expect to touch. Maybe in this case, demprep could provide an option for editing elevation cubes with spherical datums to have the right Base label value? (e.g. a new SPHERICALDATUMRADIUS parameter or similar.)

It sounds like Mars would require a more complex solution. What datum(s?) (if any) are HIRiSE DEMs usually produced relative to, as they're the ones most likely to suffer noticeably from the quantization?

[dpmayerUSGS]

@jessemapel It's not that it's been handled by GeoTiff per se, it's just that ASP is reading the projection info and assuming that the SRS describing the horizontal projection also describes a vertical datum that the pixel values are relative to.

If you wanted to make life really complicated, ISIS could try parsing compound SRSs in WKT2 format from any GDAL-supported file type ;-)

Either way, there's always going to be the risk that the geospatial header does not actually describe a vertical datum that the pixel values should be interpreted relative to. For example, I can take an HRSC DTM whose pixel values are relative to a spherical datum with radius = 3396000 meters, run map2map or gdalwarp on it and reproject it to an SRS with radius = 3396190 meters. The pixel values will still be heights relative to the 3396000 meters though unless I also do something to subtract 190 from them.

The Base/Multiplier solution doesn't completely eliminate the possibility of "de-syncing" the metadata and the pixel values as described above, but I think it makes it more difficult in practice to mess things up.

Edited to add: the Base/Multiplier solution won't work for a shape more complicated than a sphere, so parsing the A/B/C radius values in a cube or using GDAL (really, PROJ) may be necessary for a more general solution. However, I don't think WKT2 supports triaxial ellipsoids (@thareUSGS is this still true?)

[oleg-alexandrov]

Either way, there's always going to be the risk that the geospatial header does not actually describe a vertical datum that the pixel values should be interpreted relative to.

I believe it is a rather standard thing that the the projection and datum info in the geoheader of a GeoTiff are what the DEM heights in that file are relative to. Now, sure the user can choose to misuse this, but that is their problem.

[dpmayerUSGS]

@oleg-alexandrov I don't know if it's the standard thing, but it is certainly the sensible thing.

[jessemapel]

@rvwagner demprep seems like a good place to add the ability to do the base/multiplier change. It would also be helpful for non-32 bit DEMs.

[github-actions[bot]]

Thank you for your contribution!

Unfortunately, this issue hasn't received much attention lately, so it is labeled as 'stale.'

If no additional action is taken, this issue will be automatically closed in 180 days.