From 410f4bed6950d009cf75e9bd05a16741cd9adc1c Mon Sep 17 00:00:00 2001 From: "Samantha C. Niemoeller" Date: Mon, 13 Jul 2026 18:17:29 -0700 Subject: [PATCH 1/2] (not tested) layover shadow geocoding fill bugfix - convert to float32 --- .../nisar/static/layover_shadow_mask.py | 64 +++++++++++++++---- 1 file changed, 53 insertions(+), 11 deletions(-) diff --git a/python/packages/nisar/static/layover_shadow_mask.py b/python/packages/nisar/static/layover_shadow_mask.py index 4f77a54b3..d7837c4d1 100644 --- a/python/packages/nisar/static/layover_shadow_mask.py +++ b/python/packages/nisar/static/layover_shadow_mask.py @@ -195,19 +195,49 @@ def geocode_layover_shadow_mask( if geo2rdr_params is None: geo2rdr_params = {} - geocoded_layover_shadow_mask = make_scratch_gtiff( + # XXX: There's currently no way to create a "Geocode" object that operates directly + # on uint8 data with Fill value 255. 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). + # For STATIC, this shouldn't cause any loss of fidelity since the mask values + # are small and we're using nearest neighbor interpolation. + + # However, when geocoding uint8 data with `GeocodeFloat32`, the C++ `GeocodeCov` + # implementation initializes the output buffer with `nan_t_out`. Then, + # 0 * NaN evaluates to NaN for float dtypes (correct), but for integer + # dtypes it evaluates to 0 (since uint8 cannot represent NaN). This leads to + # areas outside the radar grid extent retaining this 0 fill value instead + # of the correct 255. + + # Workaround: + # 1. Convert the uint8 mask to float32, replacing 255 with NaN + # 2. Geocode the float32 raster (NaN values are preserved) + # 3. Convert back to uint8, replacing NaN with 255 + + layover_shadow_mask_float32 = make_scratch_gtiff( + shape=(radar_grid.length, radar_grid.width), + dtype=np.float32, + dir_=scratch_dir, + prefix="layover-shadow-mask-float32_", + ) + + # Convert uint8 mask to float32, replacing fill value 255 with NaN + radar_data_uint8 = np.empty((radar_grid.length, radar_grid.width), dtype=np.uint8) + layover_shadow_mask.read_all(radar_data_uint8) + radar_data_float32 = radar_data_uint8.astype(np.float32) + radar_data_float32[radar_data_uint8 == 255] = np.nan + layover_shadow_mask_float32.write_all(radar_data_float32) + + # Create output float32 raster for geocoded result + geocoded_layover_shadow_mask_float32 = make_scratch_gtiff( shape=(geo_grid.length, geo_grid.width), - dtype=np.uint8, + dtype=np.float32, dir_=scratch_dir, - prefix="geocoded-layover-shadow-mask_", + prefix="geocoded-layover-shadow-mask-float32_", ) - # 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. + # Geocode the float32 mask geocode = isce3.geocode.GeocodeFloat32() geocode.orbit = orbit geocode.ellipsoid = get_reference_ellipsoid(dem_raster) @@ -230,8 +260,8 @@ def geocode_layover_shadow_mask( geocode.geocode( radar_grid=radar_grid, - input_raster=layover_shadow_mask, - output_raster=geocoded_layover_shadow_mask, + input_raster=layover_shadow_mask_float32, + output_raster=geocoded_layover_shadow_mask_float32, dem_raster=dem_raster, output_mode=GeocodeOutputMode.INTERP, memory_mode=normalize_geocode_memory_mode(memory_mode), @@ -240,6 +270,18 @@ def geocode_layover_shadow_mask( dem_interp_method=normalize_data_interp_method(dem_interp_method), ) + # Convert geocoded float32 result back to uint8, replacing NaN with 255 + 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_", + ) + geo_data_float32 = np.empty((geo_grid.length, geo_grid.width), dtype=np.float32) + geocoded_layover_shadow_mask_float32.read_all(geo_data_float32) + geo_data_uint8 = np.where(np.isnan(geo_data_float32), 255, geo_data_float32).astype(np.uint8) + geocoded_layover_shadow_mask.write_all(geo_data_uint8) + return geocoded_layover_shadow_mask From cd41d8e46ca37ce2ee9202b83cc5965e72929d9e Mon Sep 17 00:00:00 2001 From: "Samantha C. Niemoeller" Date: Wed, 15 Jul 2026 14:43:25 -0700 Subject: [PATCH 2/2] STATIC: Use out_mask to set geocoding fill to 255 --- .../nisar/static/layover_shadow_mask.py | 86 ++++++++----------- 1 file changed, 38 insertions(+), 48 deletions(-) diff --git a/python/packages/nisar/static/layover_shadow_mask.py b/python/packages/nisar/static/layover_shadow_mask.py index d7837c4d1..f74c82cb2 100644 --- a/python/packages/nisar/static/layover_shadow_mask.py +++ b/python/packages/nisar/static/layover_shadow_mask.py @@ -196,48 +196,42 @@ def geocode_layover_shadow_mask( geo2rdr_params = {} # XXX: There's currently no way to create a "Geocode" object that operates directly - # on uint8 data with Fill value 255. 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). - # For STATIC, this shouldn't cause any loss of fidelity since the mask values - # are small and we're using nearest neighbor interpolation. - - # However, when geocoding uint8 data with `GeocodeFloat32`, the C++ `GeocodeCov` - # implementation initializes the output buffer with `nan_t_out`. Then, - # 0 * NaN evaluates to NaN for float dtypes (correct), but for integer - # dtypes it evaluates to 0 (since uint8 cannot represent NaN). This leads to - # areas outside the radar grid extent retaining this 0 fill value instead - # of the correct 255. - - # Workaround: - # 1. Convert the uint8 mask to float32, replacing 255 with NaN - # 2. Geocode the float32 raster (NaN values are preserved) - # 3. Convert back to uint8, replacing NaN with 255 - - layover_shadow_mask_float32 = make_scratch_gtiff( - shape=(radar_grid.length, radar_grid.width), - dtype=np.float32, + # 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="layover-shadow-mask-float32_", + prefix="geocoded-layover-shadow-mask_", ) - # Convert uint8 mask to float32, replacing fill value 255 with NaN - radar_data_uint8 = np.empty((radar_grid.length, radar_grid.width), dtype=np.uint8) - layover_shadow_mask.read_all(radar_data_uint8) - radar_data_float32 = radar_data_uint8.astype(np.float32) - radar_data_float32[radar_data_uint8 == 255] = np.nan - layover_shadow_mask_float32.write_all(radar_data_float32) - - # Create output float32 raster for geocoded result - geocoded_layover_shadow_mask_float32 = make_scratch_gtiff( + # 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.float32, + dtype=np.uint8, dir_=scratch_dir, - prefix="geocoded-layover-shadow-mask-float32_", + prefix="geocode-coverage-mask_", ) - # Geocode the float32 mask + # Setup the `GeocodeFloat32` object geocode = isce3.geocode.GeocodeFloat32() geocode.orbit = orbit geocode.ellipsoid = get_reference_ellipsoid(dem_raster) @@ -258,29 +252,25 @@ 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_float32, - output_raster=geocoded_layover_shadow_mask_float32, + 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), ) - # Convert geocoded float32 result back to uint8, replacing NaN with 255 - 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_", - ) - geo_data_float32 = np.empty((geo_grid.length, geo_grid.width), dtype=np.float32) - geocoded_layover_shadow_mask_float32.read_all(geo_data_float32) - geo_data_uint8 = np.where(np.isnan(geo_data_float32), 255, geo_data_float32).astype(np.uint8) - geocoded_layover_shadow_mask.write_all(geo_data_uint8) + # 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