Assign NaN to invalid pixel values during geocoding (#942)

* update the GeocodeSlc module to assign NaN to invalid pixels

* update the GeocodeInterp module to assign NaN to invalid pixels in the RTC ANF map

* set _FillValue to NaN so that NETCDF drivers assign NaN to invalid pixels

* only assign nan to _FillValue if dtype is float (_FillValue = nan) or complex (_FillValue = nan + 0j)

* substitute complex nan with nan + 1j * nan

* Update python/packages/nisar/workflows/h5_prep.py

Co-authored-by: Geoffrey M Gunter <[email protected]>

Co-authored-by: Geoffrey M Gunter <[email protected]>

Author: 2 people

cxx/isce3/geocode/GeocodeCov.cpp

@@ -760,6 +760,9 @@ inline void Geocode<T>::_interpolate(
             // (NaN, NaN)
             using T_out_real = typename isce3::real<T_out>::type;
             geoDataBlock(i, j) *= std::numeric_limits<T_out_real>::quiet_NaN();
+            if (flag_apply_rtc && out_geo_rtc != nullptr) {
+                out_geo_rtc_array(i, j) = std::numeric_limits<float>::quiet_NaN();
+            }
             continue;
         }
 
@@ -774,7 +777,6 @@ inline void Geocode<T>::_interpolate(
                     rtc_area(int(rdrY + azimuthFirstLine), int(rdrX + rangeFirstPixel));
             val /= std::sqrt(rtc_value);
             if (out_geo_rtc != nullptr) {
-#pragma omp atomic write
                 out_geo_rtc_array(i, j) = rtc_value;
             }
         }

cxx/isce3/geocode/interpolate.cpp

@@ -32,10 +32,14 @@ void isce3::geocode::interpolate(
         const double fracX = rdrX - intX;
         const double fracY = rdrY - intY;
 
-        if ((intX < chipHalf) || (intX >= (inWidth - chipHalf)))
+        if ((intX < chipHalf) || (intX >= (inWidth - chipHalf))) {
+            geoDataBlock(i, j) *= std::numeric_limits<float>::quiet_NaN();
             continue;
-        if ((intY < chipHalf) || (intY >= (inLength - chipHalf)))
+        }
+        if ((intY < chipHalf) || (intY >= (inLength - chipHalf))) {
+            geoDataBlock(i, j) *= std::numeric_limits<float>::quiet_NaN();
             continue;
+        }
 
         // Slant Range at the current output pixel
         const double rng =
@@ -46,8 +50,10 @@ void isce3::geocode::interpolate(
         const double az = radarGrid.sensingStart() +
                           radarY[i * width + j] / radarGrid.prf();
 
-        if (not dopplerLUT.contains(az, rng))
+        if (not dopplerLUT.contains(az, rng)) {
+            geoDataBlock(i, j) *= std::numeric_limits<float>::quiet_NaN();
             continue;
+        }
 
         // Evaluate Doppler at current range and azimuth time
         const double dop =

python/packages/nisar/workflows/gcov.py

@@ -536,6 +536,10 @@ def _save_hdf5_dataset(ds_filename, h5py_obj, root_path,
 
         if fill_value is not None:
             dset.attrs.create('_FillValue', data=fill_value)
+        elif 'cfloat' in gdal.GetDataTypeName(raster.datatype()).lower():
+            dset.attrs.create('_FillValue', data=np.nan + 1j * np.nan)
+        elif 'float' in gdal.GetDataTypeName(raster.datatype()).lower():
+            dset.attrs.create('_FillValue', data=np.nan)
 
         if stats_vector is not None:
             stats_obj = stats_vector[band]

python/packages/nisar/workflows/h5_prep.py

@@ -1214,7 +1214,7 @@ def _get_raster_from_hdf5_ds(group, ds_name, dtype, shape,
         del group[ds_name]
 
     # create dataset
-    dset = group.create_dataset(ds_name, dtype=np.float64, shape=shape)
+    dset = group.create_dataset(ds_name, dtype=dtype, shape=shape)
 
     if zds is not None:
         dset.dims[0].attach_scale(zds)
@@ -1239,6 +1239,10 @@ def _get_raster_from_hdf5_ds(group, ds_name, dtype, shape,
 
     if fill_value is not None:
         dset.attrs.create('_FillValue', data=fill_value)
+    elif np.issubdtype(dtype, np.floating):
+        dset.attrs.create('_FillValue', data=np.nan)
+    elif np.issubdtype(dtype, np.complexfloating):
+        dset.attrs.create('_FillValue', data=np.nan + 1j * np.nan)
 
     if valid_min is not None:
         dset.attrs.create('valid_min', data=valid_min)