Skip to content
Draft
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 38 additions & 6 deletions python/packages/nisar/static/layover_shadow_mask.py
Original file line number Diff line number Diff line change
Expand Up @@ -195,19 +195,43 @@ def geocode_layover_shadow_mask(
if geo2rdr_params is None:
geo2rdr_params = {}

# XXX: There's currently no way to create a "Geocode" object that operates directly
# on uint8 data. Instead, the way this is done in the InSAR geocoding workflow is to
# create a `GeocodeFloat32` object, resulting in conversion of the input mask values
# to float32 during geocoding (and then conversion back to uint8 when they're
# written to the output uint8 raster). For the layover/shadow mask, this shouldn't
# cause any loss of fidelity since the mask values are small and we're using
# nearest neighbor interpolation.

# However, when geocoding a `GeocodeFloat32`, the C++ `GeocodeCov` internally
# uses NaN as a fill value for the area outside the radar grid extent.
# uint8 cannot represent NaN, which leads to the geocoding fill area
# incorrectly getting the value 0. Because 0 is a valid value in the layover
# shadow mask, this is unacceptable; the geocoding fill pixels must be 255
# per the spec.

# So, we'll use geocode's `out_mask` parameter, where out_mask==255 indicates
# pixels outside radar grid extent.
# We'll use this coverage information to mask the final geocoded layover
# shadow mask, setting its geocoding fill pixels' value to 255.

# Construct the output uint8 Geotiff raster
geocoded_layover_shadow_mask = make_scratch_gtiff(
shape=(geo_grid.length, geo_grid.width),
dtype=np.uint8,
dir_=scratch_dir,
prefix="geocoded-layover-shadow-mask_",
)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
geocoded_layover_shadow_mask.fill(255)

This is what I mean. Presumably only the valid pixels get updated during processing, so if you initialize them all to 255 in the beginning then you don't have to muck around with temporary mask layers later.

I don't think we have a isce3.io.Raster.fill method yet, but it should be easy to implement. To test the approach without committing to that, you could start with

    geocoded_layover_shadow_mask.set_block(
        np.s_[:, :],
        np.full(geocoded_layover_shadow_mask.shape, 255, np.uint8))

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nemo794 Does this fix seem simpler? I can either run tests for your current commit or would you prefer to try Brian's suggestion instead?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bhawkins Thank you for clarifying! You mentioned the word "Presumably", and that unfortunately seems to be the gotcha. :/

That said, I did test your fix in hopes it would work. Unfortunately, it did not resolve the bug -- the geocoding fill area was populated with 0 values again.

I had Claude sleuth around ISCE3 to determine the reason, and here's its summary:

The `geocode.geocode()` call writes the entire `geoDataBlock` to the output raster, 
overwriting ALL pixels including the pre-filled 255s. The C++ code doesn't check 
"is this pixel already 255?" - it just blindly writes the entire block, here:

outputRaster.setBlock(geoDataBlock.data(), 0, lineStart,
geogrid.width(), geoBlockLength, band + 1);

Detailed breakdown from Claude:

The C++ Geocoding Process, from GeocodeCov.cpp lines 863-877:

// Initialize geoDataBlock with NaN (becomes 0 for uint8)
geoDataBlock.fill(nan_t_out);

_interpolate(...);  // <-- Writes ONLY to valid pixels

// Write the block to output raster
outputRaster.setBlock(geoDataBlock.data(), ...);

When you pre-fill the output raster with 255:

  • Before geocode: geocoded_layover_shadow_mask = all 255s
  • During geocode:
    • C++ creates internal geoDataBlock filled with nan_t_out (0 for uint8)
    • _interpolate() writes radar mask values (0-3) only to pixels inside radar grid
    • Pixels outside radar grid remain as 0 in geoDataBlock
  • Writing to output: C++ calls setBlock() on the output raster
    • This overwrites the pre-filled 255s with the geoDataBlock contents
    • So outside radar grid: your 255s → 0 (the bug)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A different design could maybe be to add a new fill_value parameter into GeocodeCov, so that internally it fills geoDataBlock with something other than NaN...?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. Yeah, it looks like nan_t_out ends up as zero in this case (godbolt).

Just to recap my suggestion from today's meeting, I'd do it in two parts:

  1. Add a new fill_value parameter to the function(s) that currently have a hard-coded nan_t_out, allowing a motivated user to control it.
  2. Populate the default with the values we want for NISAR static layers. That might look something like the following:
#include <cstdint>
#include <limits>
#include <type_traits>

template<typename T>
constexpr T get_default_fill_value() {
    if constexpr (std::is_same_v<T, float>) {
        return std::numeric_limits<float>::quiet_NaN();
    } else if constexpr (std::is_same_v<T, double>) {
        return std::numeric_limits<double>::quiet_NaN();
    } else if constexpr (std::is_same_v<T, unsigned char>) {
        return static_cast<unsigned char>(0xff);
    } else if constexpr (std::is_same_v<T, unsigned short>) {
        return static_cast<unsigned short>(0xffff);  // ... or 0xff maybe?
    } else {
        static_assert(!sizeof(T), "get_default_fill_value: unsupported type");
    }
}

// prototype in header
template<typename T>
void geocode(
    // ... lots of parameters
    const T fill_value = get_default_fill_value<T>());  // fn call is okay because it's constexpr

That way we get the desired behavior right away and we can work on actually percolating fill_value up through the rest of the API as time permits.

In the meeting I was thinking about using std::optional<T> fill_value = nullptr which would let you choose the default with runtime logic. But you can just do it with constexpr logic, so that's not actually necessary and you can just stick it in the header.

# XXX: There's currently no way to create a "Geocode" object that operates directly
# on uint8 data. Instead, the way this is done in the InSAR geocoding workflow is to
# create a `GeocodeFloat32` object, resulting in conversion of the input mask values
# to float32 during geocoding (and then conversion back to uint8 when they're
# written to the output raster). This shouldn't cause any loss of fidelity since the
# mask values are small and we're using nearest neighbor interpolation.
# Create output mask to indicate which geocoded pixels are outside radar grid extent
out_mask = make_scratch_gtiff(
shape=(geo_grid.length, geo_grid.width),
dtype=np.uint8,
dir_=scratch_dir,
prefix="geocode-coverage-mask_",
)

# Setup the `GeocodeFloat32` object
geocode = isce3.geocode.GeocodeFloat32()
geocode.orbit = orbit
geocode.ellipsoid = get_reference_ellipsoid(dem_raster)
Expand All @@ -228,18 +252,26 @@ def geocode_layover_shadow_mask(
if "maxiter" in geo2rdr_params:
geocode.numiter_geo2rdr = geo2rdr_params["maxiter"]

# Geocode uint8 layover shadow mask, with the out_mask
geocode.geocode(
radar_grid=radar_grid,
input_raster=layover_shadow_mask,
output_raster=geocoded_layover_shadow_mask,
dem_raster=dem_raster,
output_mode=GeocodeOutputMode.INTERP,
out_mask=out_mask,
memory_mode=normalize_geocode_memory_mode(memory_mode),
min_block_size=min_block_size,
max_block_size=max_block_size,
dem_interp_method=normalize_data_interp_method(dem_interp_method),
)

# Set geocoding fill values to 255
mask_data = geocoded_layover_shadow_mask.get_block((slice(None), slice(None)))
out_mask_data = out_mask.get_block((slice(None), slice(None)))
mask_data[out_mask_data == 255] = 255
geocoded_layover_shadow_mask.set_block((slice(None), slice(None)), mask_data)

return geocoded_layover_shadow_mask


Expand Down
Loading