Merge pull request opera-adt#110 from hb6688/dswx-ni-beta-point

Dswx ni beta point

README.md

@@ -63,7 +63,7 @@ Then, from inside the cloned repository, build the Docker image:
 Load the Docker container image onto your computer:
 
 ```bash
-docker load -i docker/dockerimg_dswx_s1_final_1.0.tar
+docker load -i docker/dockerimg_dswx_s1_final_patch_1.1.tar
 ```
 
 See DSWx-SAR Science Algorithm Software (SAS) User Guide for instructions on processing via Docker.

build_docker.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 
 IMAGE=opera/dswx-s1
-tag=final_1.0
+tag=final_patch_1.1
 echo "IMAGE is $IMAGE:$tag"
 
 # fail on any non-zero exit codes

docker/Dockerfile

@@ -1,8 +1,8 @@
 FROM oraclelinux:8
 
 LABEL author="OPERA ADT" \
-    description="DSWX-SAR final release R4" \
-    version="1.0-final"
+    description="DSWX-SAR final patch release R4.1" \
+    version="1.1-final-patch"
 
 RUN yum -y update &&\
     yum -y install curl &&\

src/dswx_sar/defaults/algorithm_parameter_s1.yaml

@@ -198,7 +198,7 @@ runconfig:
             number_cpu: 1
 
         refine_with_bimodality:
-            minimum_pixel: 4
+            minimum_pixel: 40
             lines_per_block: 500
             number_cpu: 1
             thresholds:

src/dswx_sar/defaults/dswx_s1.yaml

@@ -12,6 +12,9 @@ runconfig:
             #            for open water
             input_file_path:
 
+            # Specify the MGRS tile collection ID
+            input_mgrs_collection_id:
+
         dynamic_ancillary_file_group:
             # Digital elevation model (Required)
             dem_file:

src/dswx_sar/dswx_ni_runconfig.py

@@ -6,6 +6,7 @@
 import sys
 import logging
 from types import SimpleNamespace
+import numpy as np
 
 import yamale
 from ruamel.yaml import YAML
@@ -33,6 +34,18 @@
     'SH_POL': ['HH'],
     }
 
+# 2nd dictionary is for single frame only
+DSWX_NI_SINGLE_FRAME_POL_DICT = {
+    'SH_POL': ['HH'],
+    'SV_POL': ['VV'],
+    'DH_POL': ['HH', 'HV'],
+    'DV_POL': ['VV', 'VH'],
+    'CO_POL': ['HH', 'VV'],
+    'CROSS_POL': ['HV', 'VH'],
+    'QP_POL': ['HH', 'VV', 'HV', 'VH'],
+    'DV_SH_POL': ['VV', 'VH', 'HH'],
+    'DH_SV_POL': ['HH', 'HV', 'VV'],
+    }
 
 def _get_parser():
     parser = argparse.ArgumentParser(
@@ -190,19 +203,185 @@ def check_file_path(file_path: str) -> None:
             raise FileNotFoundError(err_str)
 
 
-def get_pol_rtc_hdf5(path):
+def get_pol_rtc_hdf5(input_rtc, freq_group):
     # convenience function to get polarization from RTC file path
     # basename separates file name from directory in path string
     # splitext removes the file extension from basename
     # split('_')[-1] gets polarization
-    path_pol = '/science/LSAR/GCOV/grids/frequencyA/listOfPolarizations'
-    with h5py.File(path) as src:
+    path_pol = f'/science/LSAR/GCOV/grids/frequency{freq_group}/listOfPolarizations'
+
+    with h5py.File(input_rtc) as src:
         pols = src[path_pol][()]
         pols = [pol.decode('utf-8') for pol in pols]
 
     return pols
 
 
+def get_freq_rtc_hdf5(input_rtc):
+    path_freq = f'/science/LSAR/identification/listOfFrequencies'
+
+    with h5py.File(input_rtc) as src_h5:
+        freq_group_list = src_h5[path_freq][()]
+        freq = [freq_group.decode('utf-8') for freq_group in freq_group_list]
+
+    return freq
+
+
+def check_rtc_frequency(input_h5_list):
+    """Read Frequency group information from input RTC file(s)
+
+    Parameters
+    ----------
+    input_h5_list : list of strings
+        input RTC files in a list
+
+    Returns
+    -------
+    freq_list: list
+        List of frequency group(s) in each of the files
+        e.g. A, B, or A and B
+    flag_freq_equal: bool
+        Flag that indicates whether the frequency groups are equal
+        among input files
+    """
+    num_input_files = len(input_h5_list)
+
+    # Handle the case when there is only one file
+    if num_input_files == 1:
+        freq_list = [get_freq_rtc_hdf5(input_h5_list[0])]
+        return True, freq_list  # If only one file, frequencies are trivially equal
+
+    freq_list = np.empty(num_input_files , dtype=object)
+    flag_pol_equal = True
+
+    for input_idx, input_h5 in enumerate(input_h5_list):
+        freq_list[input_idx] = get_freq_rtc_hdf5(input_h5)
+
+    for idx in range(num_input_files - 1):
+        if freq_list[idx] == freq_list[idx + 1]:
+            flag_freq_equal = True
+        else:
+            # Frequency groups between frames are different.
+            flag_freq_equal = False
+            break
+
+    return flag_freq_equal, freq_list
+
+
+def read_rtc_polarization(input_h5_list, freq_list):
+    num_input_files = len(input_h5_list)
+    pol_list = np.empty((num_input_files, 2) , dtype=object)
+
+    for input_idx, input_h5 in enumerate(input_h5_list):
+        extracted_strings = []
+        # Check to see if frequency group of an input file is empty
+        if freq_list[input_idx]:
+            for freq_idx, freq_group in enumerate(freq_list[input_idx]):
+                pol_list[input_idx, freq_idx] = get_pol_rtc_hdf5(input_h5, freq_group)
+
+    return pol_list
+
+
+def compare_rtc_polarization(pol_list):
+    num_input_files = len(pol_list)
+
+    pol_freq_a = pol_list[:, 0]
+    pol_freq_b = pol_list[:, 1]
+
+    pol_freq_a_first = pol_freq_a[0]
+    pol_freq_b_first = pol_freq_b[0]
+
+    flag_pol_freq_a_equal = True
+    flag_pol_freq_b_equal = True
+
+    # Compare Frequency A of all inputs
+    for pol in pol_freq_a:
+        if isinstance(pol_freq_a_first, list) and isinstance(pol, list):
+            if sorted(pol_freq_a_first) != sorted(pol):
+                flag_pol_freq_a_equal = False
+                break
+        elif pol_freq_a_first is None and pol is not None:
+            flag_pol_freq_a_equal = False
+            break
+        elif pol_freq_a_first is None and pol is None:
+            continue
+
+    # Compare Frequency B of all inputs
+    for pol in pol_freq_b:
+        if isinstance(pol_freq_b_first, list) and isinstance(pol, list):
+            if sorted(pol_freq_b_first) != sorted(pol):
+                flag_pol_freq_b_equal = False
+                break
+        elif pol_freq_b_first is None and pol is not None:
+            flag_pol_freq_b_equal = False
+            break
+        elif pol_freq_b_first is None and pol is None:
+            continue
+
+    return flag_pol_freq_a_equal, flag_pol_freq_b_equal
+
+def compare_rtc_bandwidth(bandwidth_list):
+    num_input_files = len(bandwidth_list)
+    bw_freq_a = bandwidth_list[:, 0]
+    bw_freq_b = bandwidth_list[:, 1]
+
+    if np.all(bw_freq_a == bw_freq_a[0]):
+        flag_bw_freq_a_equal = True
+        min_bw_freq_a = bw_freq_a[0]
+    else:
+        flag_bw_freq_a_equal = False
+
+        # Need to identify highest and lowest bandwidth (resolution)
+        bw_freq_a_valid = bw_freq_a > 5
+
+        if bw_freq_a_valid.size > 0:
+            min_bw_freq_a = np.min(bw_freq_a_valid) 
+        else:
+            min_bw_freq_a = bw_freq_a[0]
+
+    if np.all(bw_freq_b == bw_freq_b[0]):
+        flag_bw_freq_b_equal = True
+        min_bw_freq_b = bw_freq_b[0]
+    else:
+        flag_bw_freq_b_equal = False
+
+        # Need to identify highest and lowest bandwidth (resolution)
+        bw_freq_b_valid = bw_freq_b > 5
+
+        if bw_freq_b_valid.size > 0:
+            min_bw_freq_b = np.min(bw_freq_b_valid) 
+        else:
+            min_bw_freq_b = bw_freq_b[0]
+
+
+    return flag_bw_freq_a_equal, flag_bw_freq_b_equal, min_bw_freq_a, min_bw_freq_b
+
+def verify_nisar_mode(input_dir_list):
+    # Extract Frequency Groups of input files
+    flag_freq_equal, freq_list = check_rtc_frequency(input_dir_list)
+
+    # Extract polarizations of each frequency group of the input files
+    pol_list = read_rtc_polarization(input_dir_list, freq_list)
+
+    # Compare polariztions of frequency groups among input files
+    flag_pol_freq_a_equal, flag_pol_freq_b_equal = compare_rtc_polarization(pol_list)
+
+    # Determine NiSAR input RTC mode of operation
+    if flag_freq_equal and flag_pol_freq_a_equal and flag_pol_freq_b_equal:
+        nisar_uni_mode = True
+    else:
+        nisar_uni_mode = False
+
+    return flag_freq_equal, flag_pol_freq_a_equal, flag_pol_freq_b_equal, nisar_uni_mode
+
+def find_unique_polarizations(pol_list):
+    pol_list_flatten = pol_list.copy()
+    pol_list_flatten = np.concatenate(pol_list_flatten.ravel()).ravel()
+    pol_list_unique = set(pol_list_flatten)
+
+    return pol_list_unique
+
+
 def _find_polarization_from_data_dirs(input_h5_list):
     """
     This function walks through each directory in the given list of input directories,
@@ -224,9 +403,14 @@ def _find_polarization_from_data_dirs(input_h5_list):
         This list contains the polarization identifiers
         (like 'HH', 'VV', etc.) found in the filenames.
     """
+    num_input_rtc = len(input_h5_list)
     extracted_strings = []
-    for input_h5 in input_h5_list:
-        extracted_strings += get_pol_rtc_hdf5(input_h5)
+
+    for input_idx, input_h5 in enumerate(input_h5_list):
+        freq_strings = get_freq_rtc_hdf5(input_h5)
+
+        for freq_idx, input_freq in enumerate(freq_strings):
+            extracted_strings += get_pol_rtc_hdf5(input_h5, input_freq)
 
     # if nothing found raise error
     if not extracted_strings:
@@ -237,7 +421,7 @@ def _find_polarization_from_data_dirs(input_h5_list):
     return list(set(extracted_strings))
 
 
-def check_polarizations(pol_list, input_dir_list):
+def check_polarizations(pol_list, input_dir_list, DSWX_NI_POL_DICT):
     """
     Sort polarizations so that co-polarizations are preceded. This function
     identifies the common polarizations between the requested polarizations
@@ -250,6 +434,8 @@ def check_polarizations(pol_list, input_dir_list):
         List of polarizations.
     input_dir_list : list
         List of the input directories with RTC GeoTIFF files.
+    DSWX_NI_POL_DICT: dictionary
+        Dictionary which captures possible polarization scenarios for NISAR RTC inputs   
 
     Returns
     -------
@@ -312,6 +498,44 @@ def custom_sort(pol):
     return co_pol_list, cross_pol_list, sorted_pol_list, pol_mode
 
 
+def extract_bandwidth(input_dir_list):
+    """
+    This function extracts bandwidth of each frequency group of each input RTC.
+
+    Parameters
+    ----------
+    input_dir_list : list
+        List of the input RTC files.
+
+    Returns
+    -------
+    bandwidth_list
+        List of bandwidth from both frequency A and B of all input RTC files.
+    """
+
+    num_input_files = len(input_dir_list)
+    freq_path_list = '/science/LSAR/identification/listOfFrequencies'
+    bandwidth_list = np.zeros((num_input_files, 2), dtype=float)
+
+    for input_idx, input_rtc in enumerate(input_dir_list):
+        # Check if the file exists
+        if not os.path.exists(input_rtc):
+            raise FileNotFoundError(
+                f"The file '{input_rtc}' does not exist.")
+
+        with h5py.File(input_rtc, 'r') as src_h5:
+            freq_group_list = src_h5[freq_path_list][()]
+            freq_group_list = [freq_group.decode('utf-8') for freq_group in freq_group_list]
+
+            for freq_idx, freq_group in enumerate(freq_group_list):
+                bw_path = f'/science/LSAR/GCOV/metadata/sourceData/swaths/frequency{freq_group}/rangeBandwidth'
+                bandwidth = float(src_h5[bw_path][()])
+
+                bandwidth_list[input_idx, freq_idx] = bandwidth
+
+    return bandwidth_list
+
+
 def validate_group_dict(group_cfg: dict) -> None:
     """Check and validate runconfig entries.
     Parameters
@@ -432,8 +656,10 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str, args):
             cfg['runconfig']['groups']['input_file_group']['input_file_path']
         requested_pol = algorithm_cfg[
             'runconfig']['processing']['polarizations']
+
+        DSWX_NI_POL_DICT = DSWX_NI_SINGLE_FRAME_POL_DICT
         co_pol, cross_pol, pol_list, pol_mode = check_polarizations(
-            requested_pol, input_dir_list)
+            requested_pol, input_dir_list, DSWX_NI_POL_DICT)
 
         # update the polarizations
         algorithm_cfg['runconfig']['processing']['polarizations'] = pol_list
@@ -460,6 +686,28 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str, args):
                 f' has precedence over runconfig visualization "{debug_mode}"')
             sns.processing.debug_mode = args.debug_mode
 
+        # Determine NISAR input RTC mode of operation
+        (
+            flag_freq_equal, 
+            flag_pol_freq_a_equal, 
+            flag_pol_freq_b_equal, 
+            nisar_uni_mode
+        ) = verify_nisar_mode(input_dir_list)
+
+        # Update NiSAR processing mode
+        algorithm_cfg[
+            'runconfig']['processing']['nisar_uni_mode'] = nisar_uni_mode
+
+        # Extract bandwidth from input RTC
+        bandwidth_list = extract_bandwidth(input_dir_list)
+
+        (
+            flag_bw_freq_a_equal, 
+            flag_bw_freq_b_equal, 
+            min_bw_freq_a, 
+            min_bw_freq_b
+        ) = compare_rtc_bandwidth(bandwidth_list)
+
         log_file = sns.log_file
         if args.log_file is not None:
             logger.warning(
@@ -512,4 +760,4 @@ def as_dict(self):
     def to_yaml(self):
         self_as_dict = self.as_dict()
         yaml = YAML(typ='safe')
-        yaml.dump(self_as_dict, sys.stdout, indent=4)
+        yaml.dump(self_as_dict, sys.stdout, indent=4)