Merge pull request #91 from gshiroma/updates_2025

Updates 2025

Author: gshiroma

Docker/Dockerfile

@@ -2,7 +2,7 @@ FROM oraclelinux:8
 
 LABEL author="OPERA ADT" \
     description="RTC cal/val release R4" \
-    version="1.0.3-final"
+    version="1.0.4-final"
 
 RUN yum -y update &&\
     yum -y install curl &&\

Docker/requirements.txt

@@ -1,4 +1,4 @@
-python>=3.9,<3.10
+python
 cmake>=3.18
 eigen>=3.3
 fftw>=3.3
@@ -16,9 +16,12 @@ pyre>=1.11.2
 pytest
 cython
 ruamel.yaml
-scipy!=1.10.0
+scipy
 setuptools
 shapely
 yamale
 backoff
-isce3==0.15.0
+isce3
+libnetcdf
+libgdal-hdf5
+libgdal-netcdf

build_docker_image.sh

@@ -2,7 +2,7 @@
 
 REPO=opera
 IMAGE=rtc
-TAG=final_1.0.3
+TAG=final_1.0.4
 
 echo "IMAGE is $REPO/$IMAGE:$TAG"
 

src/rtc/defaults/rtc_s1.yaml

@@ -197,9 +197,6 @@ runconfig:
                 # OPTIONAL - Choices: "single_block", "geogrid", "geogrid_and_radargrid", and "auto" (default)
                 memory_mode:
 
-                # OPTIONAL - Processing upsampling factor applied to input geogrid
-                geogrid_upsampling: 1
-
                 # Save the incidence angle
                 save_incidence_angle: False
 

src/rtc/defaults/rtc_s1_static.yaml

@@ -198,9 +198,6 @@ runconfig:
                 # OPTIONAL - Choices: "single_block", "geogrid", "geogrid_and_radargrid", and "auto" (default)
                 memory_mode:
 
-                # OPTIONAL - Processing upsampling factor applied to input geogrid
-                geogrid_upsampling: 1
-
                 # Save the incidence angle
                 save_incidence_angle: True
 

src/rtc/h5_prep.py

@@ -191,7 +191,7 @@ def save_hdf5_file(hdf5_obj, output_hdf5_file, clip_max,
         del hdf5_obj[h5_ds]
     pol_list_s2 = np.array(pol_list, dtype='S2')
     dset = hdf5_obj.create_dataset(h5_ds, data=pol_list_s2)
-    dset.attrs['description'] = np.string_(
+    dset.attrs['description'] = np.bytes_(
                 'List of processed polarization layers')
 
     # save geogrid coordinates
@@ -294,21 +294,21 @@ def create_hdf5_file(product_id, output_hdf5_file, orbit, burst, cfg,
     '''
 
     hdf5_obj = h5py.File(output_hdf5_file, 'w')
-    hdf5_obj.attrs['Conventions'] = np.string_("CF-1.8")
-    hdf5_obj.attrs["contact"] = np.string_("[email protected]")
-    hdf5_obj.attrs["institution"] = np.string_("NASA JPL")
-    hdf5_obj.attrs["project"] = np.string_("OPERA")
-    hdf5_obj.attrs["reference_document"] = np.string_(
+    hdf5_obj.attrs['Conventions'] = np.bytes_("CF-1.8")
+    hdf5_obj.attrs["contact"] = np.bytes_("[email protected]")
+    hdf5_obj.attrs["institution"] = np.bytes_("NASA JPL")
+    hdf5_obj.attrs["project"] = np.bytes_("OPERA")
+    hdf5_obj.attrs["reference_document"] = np.bytes_(
         "Product Specification Document for the OPERA Radiometric"
         " Terrain-Corrected SAR Backscatter from Sentinel-1,"
         " JPL D-108758, Rev. Working Version 1, Aug 31, 2023")
 
     # product type
     product_type = cfg.groups.primary_executable.product_type
     if product_type == STATIC_LAYERS_PRODUCT_TYPE:
-        hdf5_obj.attrs["title"] = np.string_("OPERA RTC-S1-STATIC Product")
+        hdf5_obj.attrs["title"] = np.bytes_("OPERA RTC-S1-STATIC Product")
     else:
-        hdf5_obj.attrs["title"] = np.string_("OPERA RTC-S1 Product")
+        hdf5_obj.attrs["title"] = np.bytes_("OPERA RTC-S1 Product")
 
     populate_metadata_group(product_id, hdf5_obj, burst, cfg,
                             processing_datetime, is_mosaic)
@@ -320,23 +320,44 @@ def create_hdf5_file(product_id, output_hdf5_file, orbit, burst, cfg,
 
 
 def save_orbit(orbit, orbit_group, orbit_file_path):
+
+    # ensure that the orbit reference epoch has not fractional part
+    # otherwise, trancate it to seconds precision
+    orbit_reference_epoch = orbit.reference_epoch
+    if orbit_reference_epoch.frac != 0:
+        logger.warning('the orbit reference epoch is not an'
+                       ' integer number. Truncating it'
+                       ' to seconds precision and'
+                       ' updating the orbit ephemeris'
+                       ' accordingly.')
+
+        epoch = isce3.core.DateTime(orbit_reference_epoch.year,
+                                    orbit_reference_epoch.month,
+                                    orbit_reference_epoch.day,
+                                    orbit_reference_epoch.hour,
+                                    orbit_reference_epoch.minute,
+                                    orbit_reference_epoch.second)
+
+        orbit.update_reference_epoch(epoch)
+
     orbit.save_to_h5(orbit_group)
+
     # Add description attributes.
-    orbit_group["time"].attrs["description"] = np.string_(
+    orbit_group["time"].attrs["description"] = np.bytes_(
         "Time vector record. This"
         " record contains the time corresponding to position, velocity,"
         " acceleration records")
-    orbit_group["position"].attrs["description"] = np.string_(
+    orbit_group["position"].attrs["description"] = np.bytes_(
         "Position vector"
         " record. This record contains the platform position data with"
         " respect to WGS84 G1762 reference frame")
-    orbit_group["velocity"].attrs["description"] = np.string_(
+    orbit_group["velocity"].attrs["description"] = np.bytes_(
         "Velocity vector"
         " record. This record contains the platform velocity data with"
         " respect to WGS84 G1762 reference frame")
     orbit_group.create_dataset(
         'referenceEpoch',
-        data=np.string_(orbit.reference_epoch.isoformat()))
+        data=np.bytes_(orbit.reference_epoch.isoformat()))
 
     # Orbit source/type
     orbit_type = 'Undefined'
@@ -359,9 +380,12 @@ def save_orbit(orbit, orbit_group, orbit_file_path):
             orbit_type_list.append(orbit_type_individual)
         orbit_type = '; '.join(orbit_type_list)
 
-    d = orbit_group.require_dataset("orbitType", (), "S64",
-                                    data=np.string_(orbit_type))
-    d.attrs["description"] = np.string_(
+    if 'orbitType' in orbit_group:
+        del orbit_group['orbitType']
+    d = orbit_group.create_dataset("orbitType",
+                                   data=np.bytes_(orbit_type))
+
+    d.attrs["description"] = np.bytes_(
         "Type of orbit file used in processing")
 
 
@@ -1296,11 +1320,11 @@ def populate_metadata_group(product_id: str,
             continue
         if isinstance(data, str):
             dset = h5py_obj.create_dataset(
-                path_dataset_in_h5, data=np.string_(data))
+                path_dataset_in_h5, data=np.bytes_(data))
         else:
             dset = h5py_obj.create_dataset(path_dataset_in_h5, data=data)
 
-        dset.attrs['description'] = np.string_(description)
+        dset.attrs['description'] = np.bytes_(description)
 
 
 def save_hdf5_dataset(ds_filename, h5py_obj, root_path,
@@ -1341,7 +1365,11 @@ def save_hdf5_dataset(ds_filename, h5py_obj, root_path,
         logger.warning(f'WARNING Cannot open raster file: {ds_filename}')
         return
 
-    ds_name = layer_hdf5_dict[layer_name]
+    if isinstance(layer_name, str):
+        ds_name = layer_hdf5_dict[layer_name]
+    else:
+        ds_name = [layer_hdf5_dict[l] for l in layer_name]
+
     if long_name is not None:
         description = long_name
     else:
@@ -1377,25 +1405,27 @@ def save_hdf5_dataset(ds_filename, h5py_obj, root_path,
         dset = h5py_obj.create_dataset(h5_ds, data=data)
         dset.dims[0].attach_scale(yds)
         dset.dims[1].attach_scale(xds)
-        dset.attrs['grid_mapping'] = np.string_("projection")
+        dset.attrs['grid_mapping'] = np.bytes_("projection")
 
         if standard_name is not None:
-            dset.attrs['standard_name'] = np.string_(standard_name)
+            dset.attrs['standard_name'] = np.bytes_(standard_name)
 
         if long_name is not None:
-            dset.attrs['long_name'] = np.string_(long_name)
+            dset.attrs['long_name'] = np.bytes_(long_name)
 
-        dset.attrs['description'] = np.string_(description)
+        dset.attrs['description'] = np.bytes_(description)
 
         if units is not None:
-            dset.attrs['units'] = np.string_(units)
+            dset.attrs['units'] = np.bytes_(units)
 
         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)
+        elif 'byte' in gdal.GetDataTypeName(raster.datatype()).lower():
+            dset.attrs.create('_FillValue', data=255)
 
         if stats_vector is not None:
             stats_obj = stats_vector[band]

src/rtc/rtc_s1.py

@@ -769,22 +769,21 @@ def run_parallel(cfg: RunConfig, logfile_path, flag_logger_full_format):
 
             if flag_layover_shadow_mask_is_temporary:
                 temp_files_list.append(layover_shadow_mask_file)
-                layover_shadow_mask_file = None
             else:
                 burst_output_file_list.append(layover_shadow_mask_file)
                 logger.info(f'file saved: {layover_shadow_mask_file}')
 
-                # Take the layover shadow mask from HDF5 file if not exists
-                if save_secondary_layers_as_hdf5:
-                    layover_shadow_mask_file = (
-                        f'NETCDF:{burst_hdf5_in_output}:'
-                        f'{DATA_BASE_GROUP}/'
-                        f'{layer_hdf5_dict[LAYER_NAME_LAYOVER_SHADOW_MASK]}')
-
-                if save_mask:
-                    output_metadata_dict[
-                        LAYER_NAME_LAYOVER_SHADOW_MASK][1].append(
-                            layover_shadow_mask_file)
+            # Take the layover shadow mask from HDF5 file if not exists
+            if save_secondary_layers_as_hdf5:
+                layover_shadow_mask_file = (
+                    f'NETCDF:{burst_hdf5_in_output}:'
+                    f'{DATA_BASE_GROUP}/'
+                    f'{layer_hdf5_dict[LAYER_NAME_LAYOVER_SHADOW_MASK]}')
+
+            if save_mask:
+                output_metadata_dict[
+                    LAYER_NAME_LAYOVER_SHADOW_MASK][1].append(
+                        layover_shadow_mask_file)
 
             if not save_mask:
                 layover_shadow_mask_file = None

src/rtc/rtc_s1_single_job.py

@@ -1214,7 +1214,6 @@ def run_single_job(cfg: RunConfig):
         cfg.groups.processing.correction_lut_range_spacing_in_meters
 
     memory_mode = geocode_namespace.memory_mode
-    geogrid_upsampling = geocode_namespace.geogrid_upsampling
     shadow_dilation_size = geocode_namespace.shadow_dilation_size
     abs_cal_factor = geocode_namespace.abs_rad_cal
     clip_max = geocode_namespace.clip_max
@@ -1323,7 +1322,7 @@ def run_single_job(cfg: RunConfig):
 
     logger.info('Save layers:')
     logger.info(f'    {layer_names_dict[LAYER_NAME_LAYOVER_SHADOW_MASK]}:'
-                f' {save_rtc_anf}')
+                f' {save_mask}')
     logger.info(f'    RTC area normalization factor: {save_rtc_anf}')
     logger.info(f'    RTC area normalization factor Gamma0 to Beta0:'
                 f' {save_rtc_anf_gamma0_to_sigma0}')
@@ -1719,9 +1718,10 @@ def run_single_job(cfg: RunConfig):
             else:
                 burst_output_file_list.append(layover_shadow_mask_file)
                 logger.info(f'file saved: {layover_shadow_mask_file}')
-                if save_mask:
-                    output_metadata_dict[LAYER_NAME_LAYOVER_SHADOW_MASK][1].append(
-                        layover_shadow_mask_file)
+
+            if save_mask:
+                output_metadata_dict[LAYER_NAME_LAYOVER_SHADOW_MASK][1].append(
+                    layover_shadow_mask_file)
 
             if not save_mask:
                 layover_shadow_mask_file = None
@@ -1803,6 +1803,8 @@ def run_single_job(cfg: RunConfig):
             geo_obj.geogrid(geogrid.start_x, geogrid.start_y,
                             geogrid.spacing_x, geogrid.spacing_y,
                             geogrid.width, geogrid.length, geogrid.epsg)
+            
+            geogrid_upsampling = 1
 
             geo_obj.geocode(radar_grid=radar_grid,
                             input_raster=rdr_burst_raster,
@@ -1848,8 +1850,7 @@ def run_single_job(cfg: RunConfig):
             if product_type != STATIC_LAYERS_PRODUCT_TYPE:
                 output_imagery_list.append(geo_burst_filename)
 
-        if (flag_process and save_mask and
-                not save_secondary_layers_as_hdf5):
+        if flag_process and save_mask:
             set_mask_fill_value_and_ctable(layover_shadow_mask_file,
                                            geo_burst_filename)
 

src/rtc/runconfig.py

@@ -39,9 +39,6 @@ def load_parameters(cfg):
     if geocode_namespace.clip_min is None:
         geocode_namespace.clip_min = np.nan
 
-    if geocode_namespace.geogrid_upsampling is None:
-        geocode_namespace.geogrid_upsampling = 1.0
-
     if geocode_namespace.memory_mode == 'single_block':
         geocode_namespace.memory_mode = \
             isce3.core.GeocodeMemoryMode.SingleBlock

src/rtc/schemas/rtc_s1.yaml

@@ -223,9 +223,6 @@ geocoding_options:
     # Memory mode
     memory_mode: enum('auto', 'single_block', 'geogrid', 'geogrid_and_radargrid', required=False)
 
-    # Processing upsampling factor on top of the input geogrid
-    geogrid_upsampling: int(required=False)
-
     # Save the incidence angle
     save_incidence_angle: bool(required=False)