Add GSLC on-demand software options (#637)

* remove "unitless" units, remove non-exitent CF standard_names, remove NaN fill values, add GCOV default values for save_rtc and save_nlooks

* isolate GSLC chnages from "GCOV/GSLC On-Demand PR" #604

* fix unitest errors

* refactor code

* merge branch "develop" into update_gslc_on_demand (2)

* merge branch "develop" into update_gslc_on_demand (3)

* apply style changes

* Update GSLC RTC

* revert changes to the RTC

* remove blank spaces

* add comments to _baseband

* fix RTC.cpp format

* make differences between DEM and data interpolators more explicity

* fix doxygen

* rename GeocodeCov method data_interpolator() to dataInterpolator()

* merge develop into update_gslc_on_demand (2)

* add DEM interpolation method to geocode::loadDEM()

* Clean up white spaces in geocode unitest

* fix identation in the GeocodeCov module

* simplify variables declaration

* add comment to "magic number" 50

Author: gshiroma

cxx/isce3/geocode/GeocodeCov.cpp

@@ -45,6 +45,9 @@ class Geocode {
      * @param[out] output_raster       Output raster
      * @param[in]  dem_raster          Input DEM raster
      * @param[in]  output_mode         Geocode method
+     * @param[in]  flag_az_baseband_doppler Shift SLC azimuth spectrum to baseband
+     *  (using Doppler centroid) before interpolation
+     * @param[in]  flatten             Flatten the geocoded SLC
      * @param[in]  geogrid_upsampling  Geogrid upsampling (in each direction)
      * @param[in]  flag_upsample_radar_grid Double the radar grid sampling rate
      * @param[in]  flag_apply_rtc      Apply radiometric terrain correction (RTC)
@@ -54,7 +57,7 @@ class Geocode {
      * The value 0 indicates that the the exponent is based on the data type of
      * the input raster (1 for real and 2 for complex rasters).
      * @param[in]  rtc_min_value_db    Minimum value for the RTC area factor.
-     * Radar data with RTC area factor below this limit are ignored.
+     * Radar data with RTC area factor below this limit will be set to NaN.
      * @param[in]  rtc_geogrid_upsampling  Geogrid upsampling (in each
      * direction) used to compute the radiometric terrain correction RTC.
      * @param[in]  rtc_algorithm       RTC algorithm (RTC_BILINEAR_DISTRIBUTION
@@ -72,15 +75,20 @@ class Geocode {
      * @param[out] out_geo_rdr      Raster to which the radar-grid
      * positions (range and azimuth) of the geogrid pixels vertices will be
      * saved.
-     * @param[out] out_geo_dem     Raster to which the interpolated DEM
+     * @param[out] out_geo_dem    Raster to which the interpolated DEM
      * will be saved.
      * @param[out] out_nlooks          Raster to which the number of radar-grid
      * looks associated with the geogrid will be saved.
      * @param[out] out_geo_rtc         Output RTC area factor (in
      * geo-coordinates).
-     * @param[in]  in_rtc               Input RTC area factor (in slant-range).
-     * @param[out] out_rtc              Output RTC area factor (in slant-range).
-     * @param[in]  geocode_memory_mode  Select memory mode
+     * @param[in]  phase_screen_raster Phase screen to be removed before geocoding
+     * @param[in]  offset_az_raster    Azimuth offset raster (reference
+     * radar-grid geometry)
+     * @param[in]  offset_rg_raster    Range offset raster (reference
+     * radar-grid geometry)
+     * @param[in]  in_rtc              Input RTC area factor (in slant-range).
+     * @param[out] out_rtc             Output RTC area factor (in slant-range).
+     * @param[in]  geocode_memory_mode Select memory mode
      * @param[in]  min_block_size       Minimum block size (per thread)
      * @param[in]  max_block_size       Maximum block size (per thread)
      * @param[in]  dem_interp_method   DEM interpolation method
@@ -89,12 +97,15 @@ class Geocode {
     geocode(const isce3::product::RadarGridParameters& radar_grid,
             isce3::io::Raster& input_raster, isce3::io::Raster& output_raster,
             isce3::io::Raster& dem_raster,
-            geocodeOutputMode output_mode = geocodeOutputMode::AREA_PROJECTION,
+            geocodeOutputMode output_mode = geocodeOutputMode::INTERP,
+            bool flag_az_baseband_doppler = false,
+            bool flatten = false,
             double geogrid_upsampling = 1,
             bool flag_upsample_radar_grid = false,
             bool flag_apply_rtc = false,
-            isce3::geometry::rtcInputTerrainRadiometry input_terrain_radiometry =
-                    isce3::geometry::rtcInputTerrainRadiometry::BETA_NAUGHT,
+            isce3::geometry::rtcInputTerrainRadiometry
+                    input_terrain_radiometry = isce3::geometry::
+                            rtcInputTerrainRadiometry::BETA_NAUGHT,
             isce3::geometry::rtcOutputTerrainRadiometry
                     output_terrain_radiometry = isce3::geometry::
                             rtcOutputTerrainRadiometry::GAMMA_NAUGHT,
@@ -114,6 +125,9 @@ class Geocode {
             isce3::io::Raster* out_geo_dem = nullptr,
             isce3::io::Raster* out_geo_nlooks = nullptr,
             isce3::io::Raster* out_geo_rtc = nullptr,
+            isce3::io::Raster* phase_screen_raster = nullptr,
+            isce3::io::Raster* offset_az_raster = nullptr,
+            isce3::io::Raster* offset_rg_raster = nullptr,
             isce3::io::Raster* input_rtc = nullptr,
             isce3::io::Raster* output_rtc = nullptr,
             geocodeMemoryMode geocode_memory_mode = geocodeMemoryMode::AUTO,
@@ -130,12 +144,68 @@ class Geocode {
      * @param[in]  input_raster        Input raster
      * @param[out] output_raster       Output raster
      * @param[in]  dem_raster          Input DEM raster
+     * @param[in]  flag_az_baseband_doppler Shift SLC azimuth spectrum to baseband
+     *  (using Doppler centroid) before interpolation
+     * @param[in]  input_terrain_radiometry  Input terrain radiometry
+     * @param[in]  output_terrain_radiometry Output terrain radiometry
+     * @param[in]  rtc_min_value_db    Minimum value for the RTC area factor.
+     * Radar data with RTC area factor below this limit will be set to NaN.
+     * @param[in]  rtc_geogrid_upsampling  Geogrid upsampling (in each
+     * direction) used to compute the radiometric terrain correction RTC.
+     * @param[in]  rtc_algorithm       RTC algorithm (RTC_BILINEAR_DISTRIBUTION or
+     * RTC_AREA_PROJECTION)
+     * @param[in]  abs_cal_factor      Absolute calibration factor.
+     * @param[in]  clip_min            Clip (limit) minimum output values
+     * @param[in]  clip_max            Clip (limit) maximum output values
+     * @param[out] out_geo_rdr         Raster to which the radar-grid
+     * positions (range and azimuth) of the geogrid pixels centers will be
+     * saved.
+     * @param[out] out_geo_dem         Raster to which the interpolated DEM
+     * will be saved.
+     * @param[out] out_geo_rtc         Output RTC area factor (in
+     * geo-coordinates).
+     * @param[in]  flatten             Flatten the geocoded SLC
+     * @param[in]  phase_screen_raster Phase screen to be removed before geocoding
+     * @param[in]  offset_az_raster    Azimuth offset raster (reference
+     * radar-grid geometry)
+     * @param[in]  offset_rg_raster    Range offset raster (reference radar-grid
+     * geometry)
+     * @param[in]  in_rtc              Input RTC area factor (in slant-range).
+     * @param[out] out_rtc             Output RTC area factor (in slant-range).
+     * @param[in]  dem_interp_method   DEM interpolation method
      */
     template<class T_out>
-    void geocodeInterp(const isce3::product::RadarGridParameters& radar_grid,
-                       isce3::io::Raster& input_raster,
-                       isce3::io::Raster& output_raster,
-                       isce3::io::Raster& demRaster);
+    void geocodeInterp(
+            const isce3::product::RadarGridParameters& radar_grid,
+            isce3::io::Raster& input_raster, isce3::io::Raster& output_raster,
+            isce3::io::Raster& dem_raster,
+            bool flag_apply_rtc = false,
+            bool flag_az_baseband_doppler = false, bool flatten = false,
+            isce3::geometry::rtcInputTerrainRadiometry
+                    input_terrain_radiometry = isce3::geometry::
+                            rtcInputTerrainRadiometry::BETA_NAUGHT,
+            isce3::geometry::rtcOutputTerrainRadiometry
+                    output_terrain_radiometry = isce3::geometry::
+                            rtcOutputTerrainRadiometry::GAMMA_NAUGHT,
+            float rtc_min_value_db = std::numeric_limits<float>::quiet_NaN(),
+            double rtc_geogrid_upsampling =
+                    std::numeric_limits<double>::quiet_NaN(),
+            isce3::geometry::rtcAlgorithm rtc_algorithm =
+                    isce3::geometry::rtcAlgorithm::RTC_AREA_PROJECTION,
+            double abs_cal_factor = 1,
+            float clip_min = std::numeric_limits<float>::quiet_NaN(),
+            float clip_max = std::numeric_limits<float>::quiet_NaN(),
+            isce3::io::Raster* out_geo_rdr = nullptr,
+            isce3::io::Raster* out_geo_dem = nullptr,
+            isce3::io::Raster* out_geo_rtc = nullptr,
+            isce3::io::Raster* phase_screen_raster = nullptr,
+            isce3::io::Raster* offset_az_raster = nullptr,
+            isce3::io::Raster* offset_rg_raster = nullptr,
+            isce3::io::Raster* input_rtc = nullptr,
+            isce3::io::Raster* output_rtc = nullptr,
+            isce3::core::dataInterpMethod dem_interp_method =
+                    isce3::core::dataInterpMethod::BIQUINTIC_METHOD);
+
 
     /** Geocode using the area projection algorithm (adaptive multilooking)
      *
@@ -175,10 +245,10 @@ class Geocode {
      * geo-coordinates).
      * @param[in]  in_rtc              Input RTC area factor (in slant-range).
      * @param[out] out_rtc             Output RTC area factor (in slant-range).
-     * @param[in]  dem_interp_method   DEM interpolation method
      * @param[in]  min_block_size      Minimum block size (per thread)
      * @param[in]  max_block_size      Maximum block size (per thread)
      * @param[in]  geocode_memory_mode Select memory mode
+     * @param[in]  dem_interp_method   DEM interpolation method
      */
     template<class T_out>
     void geocodeAreaProj(
@@ -277,6 +347,8 @@ class Geocode {
 
     void doppler(isce3::core::LUT2d<double> doppler) { _doppler = doppler; }
 
+    void nativeDoppler(isce3::core::LUT2d<double> nativeDoppler) { _nativeDoppler = nativeDoppler; }
+
     void orbit(isce3::core::Orbit& orbit) { _orbit = orbit; }
 
     void ellipsoid(isce3::core::Ellipsoid& ellipsoid)
@@ -363,18 +435,87 @@ class Geocode {
                   isce3::core::ProjectionBase* _proj, int lineStart,
                   int blockLength, int blockWidth, double demMargin);
 
-    void _geo2rdr(const isce3::product::RadarGridParameters& radar_grid,
-                  double x, double y, double& azimuthTime, double& slantRange,
-                  isce3::geometry::DEMInterpolator& demInterp,
-                  isce3::core::ProjectionBase* proj);
-
+    std::string _get_nbytes_str(long nbytes);
+
+    int _geo2rdr(const isce3::product::RadarGridParameters& radar_grid,
+                 double x, double y, double& azimuthTime, double& slantRange,
+                 isce3::geometry::DEMInterpolator& demInterp,
+                 isce3::core::ProjectionBase* proj,
+                 float &dem_value);
+
+    /**
+     * @param[in] rdrDataBlock a basebanded block of data in radar coordinate
+     * @param[out] geoDataBlock a block of data in geo coordinates
+     * @param[in] radarX the radar-coordinates x-index of the pixels in geo-grid
+     * @param[in] radarY the radar-coordinates y-index of the pixels in geo-grid
+     * @param[in] radarBlockWidth width of the data block in radar coordinates
+     * @param[in] radarBlockLength length of the data block in radar coordinates
+     * @param[in] azimuthFirstLine azimuth time of the first sample
+     * @param[in] rangeFirstPixel  range of the first sample
+     * @param[in] interp interpolator object
+     * @param[in] radarGrid radar grid parameter
+     * @param[in] nativeDoppler 2D LUT Doppler of the SLC image
+     * @param[in] flatten flag to flatten the geocoded SLC
+     * @param[in] phase_screen_raster Phase screen raster
+     * @param[in] phase_screen_array  Phase screen array
+     * @param[in] abs_cal_factor      Absolute calibration factor.
+     * @param[in] clip_min            Clip (limit) minimum output values
+     * @param[in] clip_max            Clip (limit) maximum output values
+     * @param[in] flag_run_rtc        Flag to indicate if RTC is enabled
+     * @param[in] rtc_area            RTC area normalization factor array
+     * @param[out] out_geo_rtc        Output RTC area factor raster (in
+     * geo-coordinates)
+     * @param[out] out_geo_rtc        Output RTC area factor array (in
+     * geo-coordinates)
+     */
     template<class T_out>
-    void
-    _interpolate(isce3::core::Matrix<T_out>& rdrDataBlock,
-                 isce3::core::Matrix<T_out>& geoDataBlock,
-                 std::valarray<double>& radarX, std::valarray<double>& radarY,
-                 int rdrBlockWidth, int rdrBlockLength, int azimuthFirstLine,
-                 int rangeFirstPixel, isce3::core::Interpolator<T_out>* interp);
+    inline void _interpolate(const isce3::core::Matrix<T_out>& rdrDataBlock,
+                      isce3::core::Matrix<T_out>& geoDataBlock,
+                      const std::valarray<double>& radarX,
+                      const std::valarray<double>& radarY,
+                      const int radarBlockWidth, const int radarBlockLength,
+                      const int azimuthFirstLine, const int rangeFirstPixel,
+                      const isce3::core::Interpolator<T_out>* interp,
+                      const isce3::product::RadarGridParameters& radarGrid,
+                      const bool flag_az_baseband_doppler, const bool flatten,
+                      isce3::io::Raster* phase_screen_raster,
+                      isce3::core::Matrix<float>& phase_screen_array,
+                      double abs_cal_factor, float clip_min, float clip_max,
+                      bool flag_run_rtc, const isce3::core::Matrix<float>& rtc_area,
+                      isce3::io::Raster* out_geo_rtc, 
+                      isce3::core::Matrix<float>& out_geo_rtc_arraya);
+
+
+    /**
+     * param[in,out] data a matrix of data that needs to be base-banded in
+     * azimuth param[in] starting_range starting range of the data block
+     * param[in] sensing_start starting azimuth time of the data block
+     * param[in] range_pixel_spacing spacing of the slant range
+     * param[in] prf pulse repetition frequency
+     * param[in] doppler_lut 2D LUT of the image Doppler
+     */
+    template<class T2>
+    inline void _baseband(isce3::core::Matrix<T2>& data, const double starting_range,
+                   const double sensing_start, const double range_pixel_spacing,
+                   const double prf,
+                   const isce3::core::LUT2d<double>& doppler_lut);
+
+    /**
+     * param[in,out] data a matrix of data that needs to be base-banded in
+     * azimuth param[in] starting_range starting range of the data block
+     * param[in] sensing_start starting azimuth time of the data block
+     * param[in] range_pixel_spacing spacing of the slant range
+     * param[in] prf pulse repetition frequency
+     * param[in] doppler_lut 2D LUT of the image Doppler
+     */
+    template<class T2>
+    inline void _baseband(isce3::core::Matrix<std::complex<T2>>& data,
+                   const double starting_range, const double sensing_start,
+                   const double range_pixel_spacing, const double prf,
+                   const isce3::core::LUT2d<double>& doppler_lut);
+
+
+
 
     // isce3::core objects
     isce3::core::Orbit _orbit;
@@ -389,6 +530,9 @@ class Geocode {
     // radar grids parameters
     isce3::core::LUT2d<double> _doppler;
 
+    // native Doppler
+    isce3::core::LUT2d<double> _nativeDoppler;
+
     // start X position for the output geogrid
     double _geoGridStartX = std::numeric_limits<double>::quiet_NaN();
 

cxx/isce3/geocode/GeocodeCov.h

@@ -5,10 +5,11 @@ isce3::geometry::DEMInterpolator
 isce3::geocode::loadDEM(isce3::io::Raster& demRaster,
                        const isce3::product::GeoGridParameters& geoGrid,
                        int lineStart, int blockLength, int blockWidth,
-                       double demMargin)
+                       double demMargin, 
+                       isce3::core::dataInterpMethod demInterpMethod)
 {
     // DEM interpolator
-    isce3::geometry::DEMInterpolator demInterp;
+    isce3::geometry::DEMInterpolator demInterp(0, demInterpMethod);
 
     // the epsg code of the input DEM
     int epsgcode = demRaster.getEPSG();

cxx/isce3/geocode/loadDem.cpp

@@ -15,9 +15,11 @@ namespace isce3 { namespace geocode {
  * @param[in] blockLength length of the block of interest in the eocoded grid
  * @param[in] blockWidth  width of the block of interest in the eocoded grid
  * @param[in] demMargin  extra margin for the dem relative to the geocoded grid
+ * @param[in] demInterpMethod  DEM interpolation method
  */
 isce3::geometry::DEMInterpolator
 loadDEM(isce3::io::Raster& demRaster,
         const isce3::product::GeoGridParameters& geoGrid, int lineStart,
-        int blockLength, int blockWidth, double demMargin);
+        int blockLength, int blockWidth, double demMargin,
+        isce3::core::dataInterpMethod demInterpMethod = isce3::core::BICUBIC_METHOD);
 }} // namespace isce3::geocode

cxx/isce3/geocode/loadDem.h

@@ -1125,7 +1125,7 @@ void computeRtcBilinearDistribution(isce3::io::Raster& dem_raster,
     auto elapsed_time_milliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(
         std::chrono::high_resolution_clock::now() - start_time);
     float elapsed_time = ((float) elapsed_time_milliseconds.count()) / 1e3;
-    info << "elapsed time (RTC-BI) [s]: " << elapsed_time << pyre::journal::endl;    
+    info << "elapsed time (RTC-BI) [s]: " << elapsed_time << pyre::journal::endl;
 }
 
 isce3::error::ErrorCode loadDemFromProj(isce3::io::Raster& dem_raster,

cxx/isce3/geometry/RTC.cpp

@@ -145,7 +145,6 @@ void computeRtc(isce3::product::Product& product, isce3::io::Raster& dem_raster,
  * @param[in]  input_dop           Doppler LUT
  * @param[in]  dem_raster          Input DEM raster
  * @param[out] output_raster       Output raster
- * @param[in]  frequency           Product frequency
  * @param[in]  input_terrain_radiometry  Input terrain radiometry
  * @param[in]  output_terrain_radiometry Output terrain radiometry
  * @param[in]  rtc_area_mode       RTC area mode (AREA or AREA_FACTOR)

cxx/isce3/geometry/RTC.h

@@ -6,6 +6,8 @@
 #include <isce3/geometry/RTC.h>
 #include <isce3/io/Raster.h>
 
+#include <pybind11/stl.h>
+
 namespace py = pybind11;
 
 using isce3::core::parseDataInterpMethod;
@@ -24,6 +26,7 @@ void addbinding(py::class_<Geocode<T>>& pyGeocode)
     pyGeocode.def(py::init<>())
             .def_property("orbit", nullptr, &Geocode<T>::orbit)
             .def_property("doppler", nullptr, &Geocode<T>::doppler)
+            .def_property("native_doppler", nullptr, &Geocode<T>::nativeDoppler)
             .def_property("ellipsoid", nullptr, &Geocode<T>::ellipsoid)
             .def_property("threshold_geo2rdr", nullptr,
                           &Geocode<T>::thresholdGeo2rdr)
@@ -66,6 +69,8 @@ void addbinding(py::class_<Geocode<T>>& pyGeocode)
                     py::arg("output_raster"), py::arg("dem_raster"),
                     py::arg("output_mode") =
                             geocodeOutputMode::AREA_PROJECTION,
+                    py::arg("flag_az_baseband_doppler") = false,
+                    py::arg("flatten") = false,
                     py::arg("geogrid_upsampling") = 1,
                     py::arg("flag_upsample_radar_grid") = false,
                     py::arg("flag_apply_rtc") = false,
@@ -93,6 +98,9 @@ void addbinding(py::class_<Geocode<T>>& pyGeocode)
                     py::arg("out_geo_dem") = nullptr,
                     py::arg("out_geo_nlooks") = nullptr,
                     py::arg("out_geo_rtc") = nullptr,
+                    py::arg("phase_screen") = nullptr,
+                    py::arg("offset_az_raster") = nullptr,
+                    py::arg("offset_rg_raster") = nullptr,
                     py::arg("input_rtc") = nullptr,
                     py::arg("output_rtc") = nullptr,
                     py::arg("memory_mode") = geocodeMemoryMode::AUTO,
@@ -103,11 +111,9 @@ void addbinding(py::class_<Geocode<T>>& pyGeocode)
                     py::arg("dem_interp_method") = isce3::core::BIQUINTIC_METHOD);
 }
 
-
-
-void addbinding(pybind11::enum_<geocodeOutputMode> & pyGeocodeMode)
+void addbinding(pybind11::enum_<geocodeOutputMode>& pyGeocodeOutputMode)
 {
-    pyGeocodeMode.value("INTERP", geocodeOutputMode::INTERP)
+    pyGeocodeOutputMode.value("INTERP", geocodeOutputMode::INTERP)
             .value("AREA_PROJECTION", geocodeOutputMode::AREA_PROJECTION);
 };