Merge branch 'isce-framework:develop' into develop

Author: gshiroma

cxx/isce3/geometry/Topo.cpp

@@ -1,6 +1,7 @@
 #include "Topo.h"
 
 #include <algorithm>
+#include <cassert>
 #include <chrono>
 #include <cmath>
 #include <cstdio>
@@ -179,59 +180,76 @@ topo(Raster & demRaster, TopoLayers & layers)
         // Allocate vector for storing satellite position for each line
         std::vector<Vec3> satPosition(blockLength);
 
-        // For each line in block
-        double tline;
-        for (size_t blockLine = 0; blockLine < blockLength; ++blockLine) {
-
-            if (blockLine % std::max((int) (blockLength / 100), 1) == 0)
-                printf("\rTopo progress (block %d/%d): %d%%",
-                       (int) block + 1, (int) nBlocks,
-                       (int) (blockLine * 1e2 / blockLength)),
-                       fflush(stdout);
-
-            // Global line index
-            size_t line = lineStart + blockLine;
-
-            // Initialize orbital data for this azimuth line
-            Basis TCNbasis;
-            Vec3 pos, vel;
-            _initAzimuthLine(line, tline, pos, vel, TCNbasis);
-
-            satPosition[blockLine] = pos;
-
-            // Compute velocity magnitude
-            const double satVmag = vel.norm();
-
-            // For each slant range bin
-            #pragma omp parallel for reduction(+:totalconv)
-            for (size_t rbin = 0; rbin < _radarGrid.width(); ++rbin) {
-
-                // Get current slant range
-                const double rng = _radarGrid.slantRange(rbin);
-
-                // Get current Doppler value
-                const double dopfact = (0.5 * _radarGrid.wavelength()
-                                     * (_doppler.eval(tline, rng) / satVmag)) * rng;
-
-                // Store slant range bin data in Pixel
-                Pixel pixel(rng, dopfact, rbin);
-
-                // Initialize LLH to middle of input DEM and average height
-                Vec3 llh = demInterp.midLonLat();
-
-                // Perform rdr->geo iterations
-                int geostat = rdr2geo(
-                    pixel, TCNbasis, pos, vel, _ellipsoid, demInterp, llh,
-                    _radarGrid.lookSide(), _threshold, _numiter, _extraiter);
-                totalconv += geostat;
-
-                // Save data in output arrays
-                _setOutputTopoLayers(llh, layers, blockLine, pixel, pos, vel,
-                        TCNbasis, demInterp);
+        // Get the midpoint of the DEM block in LLH coordinates.
+        // This should always be safe since the call to `computeDEMBounds()`
+        // above loads a block of the DEM raster.
+        assert(demInterp.haveRaster());
+        const auto dem_midpoint = demInterp.midLonLat();
+
+        #pragma omp parallel
+        {
+            // Thread-local count of the total number of rdr2geo calls that
+            // converged successfully.
+            size_t totalconv_thread = 0;
+
+            // For each line in block
+            #pragma omp for
+            for (size_t blockLine = 0; blockLine < blockLength; ++blockLine) {
+                // Global line index
+                size_t line = lineStart + blockLine;
+
+                // Initialize orbital data for this azimuth line
+                Basis TCNbasis;
+                double tline;
+                Vec3 pos, vel;
+                _initAzimuthLine(line, tline, pos, vel, TCNbasis);
+
+                satPosition[blockLine] = pos;
+
+                // Compute velocity magnitude
+                const double satVmag = vel.norm();
+
+                // Initialize LLH to middle of input DEM block and average height.
+                auto llh = dem_midpoint;
+
+                for (size_t rbin = 0; rbin < _radarGrid.width(); ++rbin) {
+
+                    // Get current slant range
+                    const double rng = _radarGrid.slantRange(rbin);
+
+                    // Get current Doppler value
+                    const double dopfact = (0.5 * _radarGrid.wavelength()
+                                         * (_doppler.eval(tline, rng) / satVmag))
+                                         * rng;
+
+                    // Store slant range bin data in Pixel
+                    Pixel pixel(rng, dopfact, rbin);
+
+                    // Perform rdr->geo iterations
+                    int geostat = rdr2geo(
+                        pixel, TCNbasis, pos, vel, _ellipsoid, demInterp, llh,
+                        _radarGrid.lookSide(), _threshold, _numiter, _extraiter);
+                    totalconv_thread += geostat;
+
+                    // Save data in output arrays
+                    _setOutputTopoLayers(llh, layers, blockLine, pixel, pos, vel,
+                                         TCNbasis, demInterp);
+
+                    // If rdr2geo failed to converge, re-initialize the LLH estimate for
+                    // the next iteration. Otherwise, reuse the current solution as the
+                    // initial guess for the next range bin.
+                    if (geostat == 0) {
+                        llh = dem_midpoint;
+                    }
+                }
+            }
 
-            } // end OMP for loop pixels in block
-        } // end for loop lines in block
-        printf("\rTopo progress (block %d/%d): 100%%\n",
+            // Collect the total number of converged rdr2geo calls from among all
+            // threads.
+            #pragma omp atomic
+            totalconv += totalconv_thread;
+        }
+        printf("\rTopo progress (block %d/%d): Done\n",
                (int) block + 1, (int) nBlocks), fflush(stdout);
 
         // Compute layover/shadow masks for the block
@@ -447,9 +465,9 @@ computeDEMBounds(Raster & demRaster, DEMInterpolator & demInterp, size_t lineOff
 }
 
 void isce3::geometry::Topo::
-_setOutputTopoLayers(Vec3 & targetLLH, TopoLayers & layers, size_t line,
-                     Pixel & pixel, Vec3& pos, Vec3& vel, Basis & TCNbasis,
-                     DEMInterpolator & demInterp)
+_setOutputTopoLayers(const Vec3& targetLLH, TopoLayers & layers, size_t line,
+                     const Pixel& pixel, const Vec3& pos, const Vec3& vel,
+                     const Basis& TCNbasis, const DEMInterpolator& demInterp)
 {
     const double degrees = 180.0 / M_PI;
 

cxx/isce3/geometry/Topo.h

@@ -305,14 +305,14 @@ class isce3::geometry::Topo {
      * @param[in] TCNbasis basis for the line under consideration
      * @param[in] demInterp DEM interpolator object used to compute local slope
      */
-    void _setOutputTopoLayers(isce3::core::Vec3 &,
-                              TopoLayers &,
-                              size_t,
-                              isce3::core::Pixel &,
-                              isce3::core::Vec3& pos,
-                              isce3::core::Vec3& vel,
-                              isce3::core::Basis &,
-                              DEMInterpolator &);
+    void _setOutputTopoLayers(const isce3::core::Vec3& llh,
+                              TopoLayers& layers,
+                              size_t line,
+                              const isce3::core::Pixel& pixel,
+                              const isce3::core::Vec3& pos,
+                              const isce3::core::Vec3& vel,
+                              const isce3::core::Basis& TCNbasis,
+                              const DEMInterpolator& demInterp);
 
     /** Main entry point for the module; internal creation of topo rasters */
     template<typename T> void _topo(T& dem, const std::string& outdir);

python/extensions/pybind_isce3/geometry/rdr2geo.cpp

@@ -229,13 +229,13 @@ void addbinding(py::class_<Topo>& pyRdr2Geo)
             (meters or degrees)
         height_raster: isce3.io.Raster
             Output raster for height above ellipsoid (meters)
-        incidence_raster: isce3.io.Raster
+        incidence_angle_raster: isce3.io.Raster
             Output raster for incidence angle (degrees) computed from vertical
             at target
         heading_angle_raster: isce3.io.Raster
             Output raster for azimuth angle (degrees) computed anti-clockwise
             from EAST (Right hand rule)
-        local_incidence_raster: isce3.io.Raster
+        local_incidence_angle_raster: isce3.io.Raster
             Output raster for local incidence angle (degrees) at target
         local_psi_raster: isce3.io.Raster
             Output raster for local projection angle (degrees) at target

python/packages/nisar/products/__init__.py

@@ -11,3 +11,4 @@
     populate_dataset_attrs_from_spec,
 )
 from .projection import build_projection_dataset_attrs_dict
+from .utils import get_static_layers_data_access

python/packages/nisar/products/utils.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 from numpy.typing import ArrayLike
+import journal
 
 
 def to_bytes(s: str | ArrayLike) -> np.ndarray:
@@ -16,4 +17,60 @@ def to_bytes(s: str | ArrayLike) -> np.ndarray:
     numpy.ndarray
         The input string(s) converted to bytestrings with 'utf-8' encoding.
     """
-    return np.char.encode(s, encoding="utf-8")
+    return np.char.encode(s, encoding="utf-8")
+
+
+def get_static_layers_data_access(
+        static_layers_data_access_template: str | None,
+        granule_id: str | None) -> str:
+    """
+    Read the static layers data access template and replace the placeholder
+     "{granule_id}" with the granule ID, if provided.
+
+    Parameters
+    ----------
+    static_layers_data_access_template: str or None
+        Template string for static layers data access. If the string contains
+        the substring "{granule_id}", that substring will be replaced with
+        the contents of `granule_id`.
+        If set to `None` or an empty string, the function returns
+        "(NOT SPECIFIED)"
+        If the string is valid, i.e., it is not `None` or an empty string,
+        and it does not contain the substring "{granule_id}", then
+        `granule_id` is ignored.
+    granule_id: str or None
+        The granule ID, which will be used to replace the substring
+        "{granule_id}" in `static_layers_data_access_template`.
+        If the template string contains "{granule_id}", but `granule_id` is
+        `None`, empty string, or "(NOT SPECIFIED)", the function raises an
+        error.
+
+    Returns
+    -------
+    static_layers_data_access: str
+        The static layers data access string with "{granule_id}" placeholder
+        replaced. Returns "(NOT SPECIFIED)" if
+        `static_layers_data_access_template` is `None` or an empty
+        string.
+    """
+
+    if not static_layers_data_access_template:
+        return '(NOT SPECIFIED)'
+
+    static_layers_data_access = static_layers_data_access_template
+
+    if '{granule_id}' in static_layers_data_access_template:
+        if not granule_id or granule_id == '(NOT SPECIFIED)':
+            error_msg = ('The placeholder "{granule_id}" is included in'
+                         ' the static layers data access template,'
+                         ' but the `granule_id` was not provided'
+                         " or is invalid")
+            error_channel = journal.error('get_static_layers_data_access')
+            error_channel.log(error_msg)
+            raise ValueError(error_msg)
+
+        static_layers_data_access = \
+            static_layers_data_access_template.replace('{granule_id}',
+                                                       granule_id)
+
+    return static_layers_data_access

python/packages/nisar/products/writers/BaseL2WriterSingleInput.py

@@ -11,6 +11,7 @@
 from isce3.core.types import truncate_mantissa
 from isce3.geometry import get_near_and_far_range_incidence_angles
 from nisar.products.readers.orbit import load_orbit
+from nisar.products.utils import get_static_layers_data_access
 
 
 LEXICOGRAPHIC_BASE_POLS = ['HH', 'HV', 'VH', 'VV']
@@ -808,12 +809,17 @@ def populate_identification_l2_specific(self):
             'identification/platformName',
             default='(NOT SPECIFIED)')
 
-    def populate_ceos_analysis_ready_data_parameters_l2_common(self):
+        static_layers_data_access_runconfig = \
+            self.cfg['primary_executable']['static_layers_data_access']
 
-        self.copy_from_runconfig(
-            '{PRODUCT}/metadata/ceosAnalysisReadyData/staticLayersDataAccess',
-            'ceos_analysis_ready_data/static_layers_data_access',
-            default='(NOT SPECIFIED)')
+        static_layers_data_access = get_static_layers_data_access(
+            static_layers_data_access_runconfig, self.granule_id)
+
+        self.set_value(
+            'identification/staticLayersDataAccess',
+            static_layers_data_access)
+
+    def populate_ceos_analysis_ready_data_parameters_l2_common(self):
 
         ceos_ard_document_identifier = \
             ('https://ceos.org/ard/files/PFS/SAR/v1.0/CEOS-ARD_PFS'

share/nisar/defaults/gcov.yaml

@@ -134,6 +134,13 @@ runconfig:
             # the corresponding metadata field will be set to "(NOT SPECIFIED)".
             product_doi:
 
+            # Data access for the static layers product associated with the
+            # output product (URL or DOI). If not provided, the corresponding
+            # metadata field will be set to "(NOT SPECIFIED)".
+            # If this string contains the substring "{granule_id}", that substring
+            # will be replaced with the granule ID.
+            static_layers_data_access:
+
             # Composite release ID (CRID). If not provided, the corresponding
             # metadata field will be set to b'A10000'.
             composite_release_id:
@@ -547,10 +554,6 @@ runconfig:
         # to populate the output product's metadata
         ceos_analysis_ready_data:
 
-            # NOT YET IMPLEMENTED
-            # Data access for the static layers product associated with the output product (URL or DOI)
-            static_layers_data_access:
-
             # Estimated geometric accuracy (bias and standard deviation)
             # in meters in the X- and Y-directions over the coordinate
             # system defined by the product's EPSG code

share/nisar/defaults/gslc.yaml

@@ -67,6 +67,14 @@ runconfig:
             product_type: GSLC
             product_version: 0.1.0
             product_doi:
+
+            # Data access for the static layers product associated with the
+            # output product (URL or DOI). If not provided, the corresponding
+            # metadata field will be set to "(NOT SPECIFIED)".
+            # If this string contains the substring "{granule_id}", that substring
+            # will be replaced with the granule ID.
+            static_layers_data_access:
+
             composite_release_id:
             processing_type:
             processing_center:
@@ -327,10 +335,6 @@ runconfig:
         # to populate the output product's metadata
         ceos_analysis_ready_data:
 
-            # NOT YET IMPLEMENTED
-            # Data access for the static layers product associated with the output product (URL or DOI)
-            static_layers_data_access:
-
             # Estimated geometric accuracy (bias and standard deviation)
             # in meters in the X- and Y-directions over the coordinate
             # system defined by the product's EPSG code

share/nisar/schemas/gcov.yaml

@@ -62,6 +62,11 @@ runconfig:
             product_type: enum('GCOV')
             product_version: str('^\d+\.\d+\.\d+$', required=False)
             product_doi: str(required=False)
+
+            # Data access for the static layers product associated with the output
+            # product (URL or DOI)
+            static_layers_data_access: str(required=False)
+
             composite_release_id: regex(r'\w\d\d\d\d\d', name='CRID', required=False)
             processing_type: enum('PR', 'UR', 'OD', required=False)
             product_accuracy: enum('P', 'M', 'N', 'F', required=False)
@@ -537,11 +542,6 @@ geocode_options:
 
 ceos_analysis_ready_data_options:
 
-    # NOT YET IMPLEMENTED
-    # Data access for the static layers product associated with the output
-    # product (URL or DOI)
-    static_layers_data_access: str(required=False)
-
     # Estimated geometric accuracy (bias and standard deviation)
     # in meters in the X- and Y-directions over the coordinate
     # system defined by the product's EPSG code

share/nisar/schemas/gslc.yaml

@@ -63,6 +63,11 @@ runconfig:
             product_type: enum('GSLC')
             product_version: str('^\d+\.\d+\.\d+$', required=False)
             product_doi: str(required=False)
+
+            # Data access for the static layers product associated with the output
+            # product (URL or DOI)
+            static_layers_data_access: str(required=False)
+
             composite_release_id: regex(r'\w\d\d\d\d\d', name='CRID', required=False)
             processing_type: enum('PR', 'UR', 'OD', required=False)
             product_accuracy: enum('P', 'M', 'N', 'F', required=False)
@@ -302,10 +307,6 @@ geocode_options:
 
 ceos_analysis_ready_data_options:
 
-    # NOT YET IMPLEMENTED
-    # Data access for the static layers product associated with the output product (URL or DOI)
-    static_layers_data_access: str(required=False)
-
     # Estimated geometric accuracy (bias and standard deviation)
     # in meters in the X- and Y-directions over the coordinate
     # system defined by the product's EPSG code