Bugfix: STATIC Layover Shadow Mask geocoding fill not set to 255

[nemo794]

WIP - unit tests are in-progress. Code is ready for review/discussion.

Summary

When the STATIC workflow was first implemented in #144 , each raster layer provided complete coverage of the given track frame's extent; there was never any geocoding fill.

Then, #178 and #205 introduced the option into the STATIC workflow to only compute the science values for a smaller swath within that track frame, to e.g. reduce memory usage. This meant that several STATIC layer rasters would now include geocoding fill regions around the corners/edges of the rasters.

For raster layers that are float-valued, the geocoding fill pixels was correctly set to NaN.

However, for the layover shadow mask, the geocoding fill pixels were incorrectly set to 0. This bug is due to a limitation of GeocodeCov not being able to natively handle geocoding uint8 layers (Full details noted in the in-line code comments).

This PR resolves that issue so that the layover shadow mask layer has the geocoding fill correctly set to 255.

Technical Design

Considerations

Related Ideas to keep in mind when selecting the correct implementation design:

1. This issue also impacts GUNW's Connected Components layer, and potentially others.

• For example, if nisarqadecides to generate PNGs for connected components, it would need to implement a similar fix.
• Whatever approach is decided here will need to be adopted by other workflows.

2. Choice of Design Pattern

• Design Option 1 - use out_mask
  • Currently, this PR implements a fix by having the calling function setup an out_mask parameter, which GeocodeCov then populates with the geocoding fill area. The calling function is then responsible for applying that mask to set the geocoding fill pixel values correctly to 255.
    • Pros: reduces number of rasters required to be in memory
    • Cons: requires an extra scratch Geotiff file
• Design Option 2 - convert to float32
  • A different approach could be to 1) Convert the uint8 mask to float32, replacing 255 with NaN, then 2) Geocode the float32 raster (NaN values are preserved), then 3) Convert back to uint8, replacing NaN with 255.
    • Cons: requires add'l copies of the raster to be in memory
• Design Option 3 - wrap uint8 values
  • A trick with uint8 is that if you add 255, then instead of overflow, it will wrap values back around to e.g. 0. This works with numpy, would need more research re: robustness in other situations.
  • Pseudo code: 1) Take the uint8 mask layer and add 1 to all values. 2) Geocode; the fill values will be set to 0. 3) Add 255 to the geocoded mask. This will set geocoding fill to 255 (correct), and all other values will get wrapped back to their original correct values.
• Design Option 4 - Choose 1, 2, or 3, but implement in a light wrapper around GeocodeCov
  • Require 1-line change to all calling functions
  • Pros:
    • Does not touch GeocodeCov, which can be finnicky
    • Reasonable to maintain long-term
  • Users with uint8 datasets will need to remember to call the wrapper function
• Design Option 5 - Choose 1, 2, or 3, but implement internal to GeocodeCov
  • Should automatically resolve the issue across all products with minimal changes to their implementation
    • Pros: Reduces technical debt. Will improve user experience in the future.
    • Cons: Potentially messier refactor, which touches a core part of the ISCE3 code.

Design Implemented

This PR implements Design Option 1 for only the STATIC Layers workflow.

Open to discussion about whether a different design would be preferred.

Test Results

Tested locally. Once the design pattern is approved, @jplzhan could you please test this branch in the production pipeline?

Before this PR

After this PR

(Fill Value of 255 is a lighter grey; it is not noted in the colorbar. I'll update the plot if/when there's time)

PDF illustrating all layers, and how the radar grid varies between layers:
NISAR_L2_STATIC_023_D_067_0800_0800_19991231T235959_A10000_J_001_report.pdf

cc: @hfattahi

[nemo794]

Tagging @xhuang-jpl and @seongsujeong , because this design pattern should be implemented for GUNW's connected components layer, too.

[bhawkins]

Why not just set all pixels in the raster to 255 immediately after creation? The value 255 can be represented exactly in float32, so I don't see why conversions to/from uint8 are an issue.

[nemo794]

Why not just set all pixels in the raster to 255 immediately after creation? The value 255 can be represented exactly in float32, so I don't see why conversions to/from uint8 are an issue.

I'm not following. 🤔 "set all pixels in the raster to 255 immediately after creation" -> Which raster are you referring to? At what step in the processing?

[bhawkins]

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))

[jplzhan]

@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?

[nemo794]

@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:

isce3/cxx/isce3/geocode/GeocodeCov.cpp

Lines 894 to 895 in bdf1f6f

	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)

[nemo794]

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...?

[bhawkins]

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.

[jplzhan]

I ran tests for the current commit as is (cd41d8e) and placed the results here:

aws s3 ls s3://nisar-adt/zhan/layover_shadowmask_fix/
                           PRE 013_071/
                           PRE 015_056/
                           PRE 089_050/
                           PRE 107_123/