Merge branch 'isce-framework:develop' into develop

gshiroma · web-flow · commit a75e7f1b71f0 · 2025-08-11T12:47:42.000-07:00
diff --git a/python/packages/isce3/core/__init__.py b/python/packages/isce3/core/__init__.py
@@ -16,5 +16,6 @@
 from .poly2d import fit_bivariate_polynomial
 from . import rdr_geo_block_generator
 from .block_param_generator import BlockParam
+from .projections import is_utm
 from .serialization import load_orbit_from_h5_group
 from . import types
diff --git a/python/packages/isce3/core/projections.py b/python/packages/isce3/core/projections.py
@@ -0,0 +1,16 @@
+def is_utm(epsg: int) -> bool:
+    """
+    Check if the input EPSG code represents a UTM zone.
+
+    Parameters
+    ----------
+    epsg : int
+        The input EPSG code.
+
+    Returns
+    -------
+    bool
+        True if the coordinate reference system represented by `epsg` is a Universal
+        Transverse Mercator (UTM) projection; otherwise False.
+    """
+    return (32600 < epsg <= 32660) or (32700 < epsg <= 32760)
diff --git a/python/packages/nisar/products/__init__.py b/python/packages/nisar/products/__init__.py
@@ -9,3 +9,4 @@
     populate_global_attrs_from_spec,
     populate_dataset_attrs_from_spec,
 )
+from .projection import build_projection_dataset_attrs_dict
diff --git a/python/packages/nisar/products/insar/GOFF_writer.py b/python/packages/nisar/products/insar/GOFF_writer.py
@@ -69,7 +69,7 @@ def add_parameters_to_procinfo_group(self):
         for rslc_name in ['reference', 'secondary']:
             rslc = self[self.group_paths.ParametersPath][rslc_name]
             rslc['referenceTerrainHeight'].attrs['description'] = \
-                np.bytes_("Reference Terrain Height as a function of"
+                np.bytes_("Reference terrain height as a function of"
                            f" map coordinates for {rslc_name} RSLC")
             rslc['referenceTerrainHeight'].attrs['units'] = \
                 Units.meter
@@ -124,7 +124,7 @@ def add_grids_to_hdf5(self):
                  " where 1 is water and 0 is non-water;"
                  " the second digit represents the subswath number of that pixel in the reference RSLC;"
                  " the least-significant digit represents the subswath number of that pixel in the secondary RSLC."
-                 " A value of '0' in either subswath digit indicates an invalid sample in the corresponding RSLC"),
+                 " A value of 0 in either subswath digit indicates an invalid sample in the corresponding RSLC"),
                 grid_mapping=grids_val,
                 xds=xds,
                 yds=yds,
diff --git a/python/packages/nisar/products/insar/GUNW_writer.py b/python/packages/nisar/products/insar/GUNW_writer.py
@@ -186,7 +186,7 @@ def add_parameters_to_procinfo_group(self):
         for rslc_name in ['reference', 'secondary']:
             rslc = self[self.group_paths.ParametersPath][rslc_name]
             rslc['referenceTerrainHeight'].attrs['description'] = \
-                np.bytes_("Reference Terrain Height as a function of"
+                np.bytes_("Reference terrain height as a function of"
                            f" map coordinates for {rslc_name} RSLC")
             rslc['referenceTerrainHeight'].attrs['units'] = \
                 Units.meter
@@ -273,7 +273,7 @@ def add_grids_to_hdf5(self):
                      " where 1 is water and 0 is non-water;"
                      " the second digit represents the subswath number of that pixel in the reference RSLC;"
                      " the least-significant digit represents the subswath number of that pixel in the secondary RSLC."
-                     " A value of '0' in either subswath digit indicates an invalid sample in the corresponding RSLC"),
+                     " A value of 0 in either subswath digit indicates an invalid sample in the corresponding RSLC"),
                     grid_mapping=grids_val,
                     xds=xds,
                     yds=yds,
diff --git a/python/packages/nisar/products/insar/InSAR_L1_writer.py b/python/packages/nisar/products/insar/InSAR_L1_writer.py
@@ -502,7 +502,7 @@ def add_pixel_offsets_to_swaths_group(self):
                                                  " subswath number of that pixel in the secondary RSLC,"
                                                  " and the most significant digit represents"
                                                  " the subswath number of that pixel in the reference RSLC."
-                                                 " A value of '0' in either digit indicates an invalid sample"
+                                                 " A value of 0 in either digit indicates an invalid sample"
                                                  " in the corresponding RSLC"),
                                     fill_value=255)
             offset_group['mask'].attrs['long_name'] = np.bytes_("Valid samples subswath mask")
diff --git a/python/packages/nisar/products/insar/InSAR_base_writer.py b/python/packages/nisar/products/insar/InSAR_base_writer.py
@@ -4,6 +4,7 @@
 from typing import Any, Optional, Union
 
 import h5py
+import journal
 import numpy as np
 from isce3.core import crop_external_orbit
 from isce3.core.types import complex32, to_complex32
@@ -384,7 +385,7 @@ def add_RSLC_to_procinfo_params_group(self, rslc_name: str):
 
         reference_terrain_height = "referenceTerrainHeight"
         reference_terrain_height_description = \
-            f"Reference Terrain Height as a function of time for {rslc_name} RSLC"
+            f"Reference terrain height as a function of time for {rslc_name} RSLC"
         if reference_terrain_height in src_param_group:
             src_param_group.copy(reference_terrain_height, dst_param_group)
             dst_param_group[reference_terrain_height].attrs['description'] = \
@@ -914,6 +915,7 @@ def add_identification_to_hdf5(self):
         """
         Add the identification group to the product
         """
+        warning_channel = journal.warning('InSAR_base_writer.add_identification_to_hdf5')
         radar_band_name = self._get_band_name()
         primary_exec_cfg = self.cfg["primary_executable"]
 
@@ -930,10 +932,13 @@ def add_identification_to_hdf5(self):
             processing_type = np.bytes_('Nominal')
         elif processing_type == 'UR':
             processing_type = np.bytes_('Urgent')
-        elif processing_type == 'OD':
-            processing_type = np.bytes_('Custom')
         else:
-            processing_type = np.bytes_('Undefined')
+            processing_type = np.bytes_('Custom')
+            if processing_type != 'OD':
+                warning_channel.log(
+                    'The processing type in the runconfig is set to'
+                    f' "{processing_type}", which is not a valid value'
+                    ' for the output product metadata. Defaulting to "Custom"')
 
         # Adopt same logic as RSLC, GSLC, GCOV
         # If no condition is met, assign string from runconfig
@@ -1162,7 +1167,10 @@ def add_identification_to_hdf5(self):
             DatasetParams(
                 "processingType",
                 processing_type,
-                "Nominal (or) Urgent (or) Custom (or) Undefined",
+                'Processing pipeline used to generate this granule. ' \
+                '"Nominal": standard production system; "Urgent": time-sensitive ' \
+                'processing in response to urgent response events; ' \
+                '"Custom": user-initiated processing outside the nominal production system',
             ),
             DatasetParams(
                 "radarBand", radar_band_name,
diff --git a/python/packages/nisar/products/projection.py b/python/packages/nisar/products/projection.py
@@ -0,0 +1,97 @@
+import numpy as np
+from osgeo import osr
+from numpy.typing import ArrayLike
+
+import isce3
+from isce3.product.cf_conventions import get_grid_mapping_name
+
+
+def build_projection_dataset_attrs_dict(epsg: int) -> dict[str, ArrayLike]:
+    """
+    Get attributes describing the spatial reference system of NISAR L2 products.
+
+    Returns a dict that may be used to populate the Attributes of `projection` Datasets
+    in NISAR Level 2 products in a way that's compliant with Climate and Forecast (CF)
+    Metadata Conventions.
+
+    Parameters
+    ----------
+    epsg : int
+        The EPSG code associated with the spatial reference system.
+
+    Returns
+    -------
+    dict
+        A dict containing attributes describing the spatial reference system.
+        String-valued attributes are represented by bytestrings with 'utf-8' encoding.
+    """
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(epsg)
+
+    # Common attributes for all spatial reference systems.
+    # FIXME: Should `spatial_ref` be replaced with `crs_wkt` for CF-1.7 compliance? See
+    # https://github-fn.jpl.nasa.gov/NISAR-ADT/NISAR_PIX/issues/318.
+    attrs = dict(
+        epsg_code=epsg,
+        spatial_ref=np.bytes_(srs.ExportToWkt()),
+        grid_mapping_name=np.bytes_(get_grid_mapping_name(srs)),
+        semi_major_axis=6378137.0,
+        inverse_flattening=298.257223563,
+        ellipsoid=np.bytes_("WGS84"),
+    )
+
+    if epsg == 4326:  # WGS 84
+        attrs["longitude_of_prime_meridian"] = 0.0
+        return attrs
+
+    # Common attributes for all projected spatial reference systems.
+    attrs["false_easting"] = srs.GetProjParm(osr.SRS_PP_FALSE_EASTING)
+    attrs["false_northing"] = srs.GetProjParm(osr.SRS_PP_FALSE_NORTHING)
+    attrs["longitude_of_projection_origin"] = srs.GetProjParm(
+        osr.SRS_PP_LONGITUDE_OF_ORIGIN
+    )
+
+    if epsg == 3413:  # Polar Stereographic (North)
+        attrs["latitude_of_projection_origin"] = 90.0
+        attrs["standard_parallel"] = 70.0
+        attrs["straight_vertical_longitude_from_pole"] = -45.0
+        return attrs
+
+    if epsg == 3031:  # Polar Stereographic (South)
+        attrs["latitude_of_projection_origin"] = -90.0
+        attrs["standard_parallel"] = -71.0
+        attrs["straight_vertical_longitude_from_pole"] = 0.0
+        return attrs
+
+    # Attribute for non-Polar-Stereographic spatial reference systems.
+    attrs["latitude_of_projection_origin"] = srs.GetProjParm(
+        osr.SRS_PP_LATITUDE_OF_ORIGIN
+    )
+
+    if isce3.core.is_utm(epsg):
+        attrs["utm_zone_number"] = epsg % 100
+        attrs["longitude_of_central_meridian"] = srs.GetProjParm(
+            osr.SRS_PP_CENTRAL_MERIDIAN
+        )
+        attrs["scale_factor_at_central_meridian"] = srs.GetProjParm(
+            osr.SRS_PP_SCALE_FACTOR
+        )
+        return attrs
+
+    if epsg == 6933:  # EASE-Grid 2.0
+        attrs["longitude_of_central_meridian"] = 0.0
+        attrs["standard_parallel"] = 30.0
+        return attrs
+
+    # FIXME: This CRS uses the GRS80 reference ellipsoid -- not WGS 84. The ellipsoid
+    # parameters defined above are not valid for this EPSG code. LAEA Europe is not
+    # supported anywhere else in ISCE3. Should we remove this? See
+    # https://github.com/isce-framework/isce3/issues/72.
+    if epsg == 3035:  # LAEA Europe
+        attrs["standard_parallel"] = -71.0
+        attrs["straight_vertical_longitude_from_pole"] = 0.0
+        return attrs
+
+    raise NotImplementedError(
+        f"EPSG {epsg} waiting for implementation / not supported in ISCE3"
+    )
diff --git a/python/packages/nisar/workflows/h5_prep.py b/python/packages/nisar/workflows/h5_prep.py
@@ -7,14 +7,13 @@
 
 import h5py
 import isce3
-from isce3.product.cf_conventions import get_grid_mapping_name
 import journal
 import numpy as np
 from isce3.core.types import complex32, to_complex32
 from isce3.io import HDF5OptimizedReader, optimize_chunk_size
 from nisar.h5 import cp_h5_meta_data
+from nisar.products.projection import build_projection_dataset_attrs_dict
 from nisar.products.readers import SLC
-from osgeo import osr
 
 
 def get_dataset_output_options(cfg: dict):
@@ -769,100 +768,19 @@ def set_get_geo_info(hdf5_obj, root_ds, geo_grid, z_vect=None,
                                       shape=(),
                                       dtype=np.uint32,
                                       data=epsg_code)
-    # Set up osr for wkt
-    srs = osr.SpatialReference()
-    srs.ImportFromEPSG(epsg_code)
+
+    projds_attrs = build_projection_dataset_attrs_dict(epsg_code)
+    projds.attrs.update(projds_attrs)
 
     # Add projection description
     projds.attrs['description'] = np.bytes_('Product map grid projection: EPSG code, '
                                              'with additional projection information as HDF5 Attributes')
 
-    # WGS84 ellipsoid
-    projds.attrs['semi_major_axis'] = 6378137.0
-    projds.attrs['inverse_flattening'] = 298.257223563
-    projds.attrs['ellipsoid'] = np.bytes_("WGS84")
-
-    # Additional fields
-    projds.attrs['epsg_code'] = epsg_code
-
-    # CF 1.7+ requires this attribute to be named "crs_wkt"
-    # spatial_ref is old GDAL way. Using that for testing only.
-    # For NISAR replace with "crs_wkt"
-    projds.attrs['spatial_ref'] = np.bytes_(srs.ExportToWkt())
-
-    # Here we have handcoded the attributes for the different cases
-    # Recommended method is to use pyproj.CRS.to_cf() as shown above
-    # To get complete set of attributes.
-
-    sr = osr.SpatialReference()
-    sr.ImportFromEPSG(epsg_code)
-
-    projds.attrs['grid_mapping_name'] = np.bytes_(get_grid_mapping_name(sr))
-
-    # Set up units
-    # Geodetic latitude / longitude
-    if epsg_code == 4326:
-        # Set up grid mapping
-        projds.attrs['longitude_of_prime_meridian'] = 0.0
-    else:
-        # UTM zones
-        if ((epsg_code > 32600 and
-             epsg_code < 32661) or
-                (epsg_code > 32700 and
-                 epsg_code < 32761)):
-            # Set up grid mapping
-            projds.attrs['utm_zone_number'] = epsg_code % 100
-            projds.attrs["longitude_of_central_meridian"] = srs.GetProjParm(
-                osr.SRS_PP_CENTRAL_MERIDIAN)
-            projds.attrs["scale_factor_at_central_meridian"] = srs.GetProjParm(
-                osr.SRS_PP_SCALE_FACTOR)
-
-        # Polar Stereo North
-        elif epsg_code == 3413:
-            # Set up grid mapping
-            projds.attrs['latitude_of_projection_origin'] = 90.0
-            projds.attrs['standard_parallel'] = 70.0
-            projds.attrs['straight_vertical_longitude_from_pole'] = -45.0
-
-        # Polar Stereo south
-        elif epsg_code == 3031:
-            # Set up grid mapping
-            projds.attrs['latitude_of_projection_origin'] = -90.0
-            projds.attrs['standard_parallel'] = -71.0
-            projds.attrs['straight_vertical_longitude_from_pole'] = 0.0
-
-        # EASE 2 for soil moisture L3
-        elif epsg_code == 6933:
-            # Set up grid mapping
-            projds.attrs['longitude_of_central_meridian'] = 0.0
-            projds.attrs['standard_parallel'] = 30.0
-
-        # Europe Equal Area for Deformation map (to be implemented in isce3)
-        elif epsg_code == 3035:
-            # Set up grid mapping
-            projds.attrs['standard_parallel'] = -71.0
-            projds.attrs['straight_vertical_longitude_from_pole'] = 0.0
-
-        else:
-            raise NotImplementedError(
-                f'EPSG {epsg_code} waiting for implementation / not supported in ISCE3')
-
+    if epsg_code != 4326:
         # Setup common parameters
         xds.attrs['long_name'] = np.bytes_("X coordinate of projection")
         yds.attrs['long_name'] = np.bytes_("Y coordinate of projection")
 
-        projds.attrs['false_easting'] = sr.GetProjParm(osr.SRS_PP_FALSE_EASTING)
-        projds.attrs['false_northing'] = sr.GetProjParm(
-            osr.SRS_PP_FALSE_NORTHING)
-
-        projds.attrs['longitude_of_projection_origin'] = sr.GetProjParm(
-            osr.SRS_PP_LONGITUDE_OF_ORIGIN)
-
-        if epsg_code not in [3413, 3031]:
-            projds.attrs['latitude_of_projection_origin'] = sr.GetProjParm(
-                osr.SRS_PP_LATITUDE_OF_ORIGIN)
-
-
     if z_vect is not None:
         return zds, yds, xds
     return yds, xds
@@ -924,25 +842,33 @@ def add_radar_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, geogrid,
         cube_group, 'losUnitVectorX', np.float32, cube_shape,
         zds=zds, yds=yds, xds=xds,
         long_name='LOS unit vector X',
-        descr='East component of unit vector of LOS from target to sensor',
+        descr='East component of the line-of-sight (LOS) unit vector, defined from ' \
+              'the target to the sensor, expressed in the east-north-up (ENU) coordinate ' \
+              'system with its origin at the target location',
         units='1', valid_min=-1.0, valid_max=1.0, **create_dataset_kwargs)
     los_unit_vector_y_raster = _get_raster_from_hdf5_ds(
         cube_group, 'losUnitVectorY', np.float32, cube_shape,
         zds=zds, yds=yds, xds=xds,
         long_name='LOS unit vector Y',
-        descr='North component of unit vector of LOS from target to sensor',
+        descr='North component of the line-of-sight (LOS) unit vector, defined from ' \
+              'the target to the sensor, expressed in the east-north-up (ENU) coordinate ' \
+              'system with its origin at the target location',
         units='1', valid_min=-1.0, valid_max=1.0, **create_dataset_kwargs)
     along_track_unit_vector_x_raster = _get_raster_from_hdf5_ds(
         cube_group, 'alongTrackUnitVectorX', np.float32, cube_shape,
         zds=zds, yds=yds, xds=xds,
         long_name='Along-track unit vector X',
-        descr='East component of unit vector along ground track',
+        descr='East component of the along-track unit vector at the target location, ' \
+              'expressed in the east-north-up (ENU) coordinate system and projected ' \
+              'onto the horizontal plane (i.e., excluding the up component)',
         units='1', valid_min=-1.0, valid_max=1.0, **create_dataset_kwargs)
     along_track_unit_vector_y_raster = _get_raster_from_hdf5_ds(
         cube_group, 'alongTrackUnitVectorY', np.float32, cube_shape,
         zds=zds, yds=yds, xds=xds,
         long_name='Along-track unit vector Y',
-        descr='North component of unit vector along ground track',
+        descr='North component of the along-track unit vector at the target location, ' \
+              'expressed in the east-north-up (ENU) coordinate system and projected ' \
+              'onto the horizontal plane (i.e., excluding the up component)',
         units='1', valid_min=-1.0, valid_max=1.0, **create_dataset_kwargs)
     elevation_angle_raster = _get_raster_from_hdf5_ds(
         cube_group, 'elevationAngle', np.float32, cube_shape,
diff --git a/tests/python/packages/CMakeLists.txt b/tests/python/packages/CMakeLists.txt
@@ -46,6 +46,7 @@ nisar/products/readers/raw.py
 nisar/products/readers/rslc.py
 nisar/products/readers/rslc_cal.py
 nisar/products/product_spec.py
+nisar/products/projection.py
 nisar/workflows/crossmul.py
 nisar/workflows/estimate_abscal_factor.py
 nisar/workflows/faraday_rot_angle_from_rslc.py
diff --git a/tests/python/packages/nisar/products/projection.py b/tests/python/packages/nisar/products/projection.py

Original file line number	Diff line number	Diff line change
`@@ -9,3 +9,4 @@`
`9`	`9`	`populate_global_attrs_from_spec,`
`10`	`10`	`populate_dataset_attrs_from_spec,`
`11`	`11`	`)`
	`12`	`+from .projection import build_projection_dataset_attrs_dict`