Merge pull request #74 from SAR-ARD/feature/sigma0

sigma0 processing and annotation layer configuration
SAR-ARD · Dec 29, 2022 · 2d48c8f · 2d48c8f
2 parents 67ab40d + f43228d
commit 2d48c8f
Show file tree

Hide file tree

Showing 12 changed files with 465 additions and 145 deletions.
diff --git a/S1_NRB/config.py b/S1_NRB/config.py
@@ -22,16 +22,17 @@ def get_config(config_file, proc_section='PROCESSING'):
     if not os.path.isfile(config_file):
         raise FileNotFoundError("Config file {} does not exist.".format(config_file))
 
-    parser = configparser.ConfigParser(allow_no_value=True, converters={'_datetime': _parse_datetime,
+    parser = configparser.ConfigParser(allow_no_value=True, converters={'_annotation': _parse_annotation,
+                                                                        '_datetime': _parse_datetime,
                                                                         '_tile_list': _parse_tile_list})
     parser.read(config_file)
     out_dict = {}
 
     # PROCESSING section
     allowed_keys = ['mode', 'aoi_tiles', 'aoi_geometry', 'mindate', 'maxdate', 'acq_mode',
-                    'work_dir', 'scene_dir', 'rtc_dir', 'tmp_dir', 'wbm_dir',
+                    'work_dir', 'scene_dir', 'rtc_dir', 'tmp_dir', 'wbm_dir', 'measurement',
                     'db_file', 'kml_file', 'dem_type', 'gdal_threads', 'log_dir', 'nrb_dir',
-                    'etad', 'etad_dir', 'product']
+                    'etad', 'etad_dir', 'product', 'annotation']
     try:
         proc_sec = parser[proc_section]
     except KeyError:
@@ -90,8 +91,22 @@ def get_config(config_file, proc_section='PROCESSING'):
         if k == 'product':
             allowed = ['GRD', 'SLC']
             assert v in allowed, "Parameter '{}': expected to be one of {}; got '{}' instead".format(k, allowed, v)
+        if k == 'measurement':
+            allowed = ['gamma', 'sigma']
+            assert v in allowed, "Parameter '{}': expected to be one of {}; got '{}' instead".format(k, allowed, v)
+        if k == 'annotation':
+            v = proc_sec.get_annotation(k)
         out_dict[k] = v
 
+    # use previous defaults for measurement and annotation if they have not been defined
+    if 'measurement' not in out_dict.keys():
+        out_dict['measurement'] = 'gamma'
+    if 'annotation' not in out_dict.keys():
+        if out_dict['measurement'] == 'gamma':
+            out_dict['annotation'] = ['dm', 'ei', 'id', 'lc', 'li', 'np', 'gs']
+        else:
+            out_dict['annotation'] = ['dm', 'ei', 'id', 'lc', 'li', 'np', 'sg']
+
     assert any([out_dict[k] is not None for k in ['aoi_tiles', 'aoi_geometry']])
 
     # METADATA section
@@ -110,6 +125,21 @@ def get_config(config_file, proc_section='PROCESSING'):
     return out_dict
 
 
+def _parse_annotation(s):
+    """Custom converter for configparser:
+    https://docs.python.org/3/library/configparser.html#customizing-parser-behaviour"""
+    if s in ['', 'None']:
+        return None
+    annotation_list = s.replace(' ', '').split(',')
+    allowed = ['dm', 'ei', 'em', 'id', 'lc', 'li', 'np', 'gs', 'sg']
+    for annotation in annotation_list:
+        if annotation not in allowed:
+            msg = "Parameter 'annotation': Error while parsing to list; " \
+                  "annotation '{}' is not supported. Allowed keys:\n{}"
+            raise ValueError(msg.format(annotation, allowed))
+    return annotation_list
+
+
 def _parse_datetime(s):
     """Custom converter for configparser:
     https://docs.python.org/3/library/configparser.html#customizing-parser-behaviour"""
@@ -171,8 +201,6 @@ def snap_conf(config):
                         'SM': 10,
                         'EW': 20}[config['acq_mode']],
             'allow_res_osv': True,
-            'export_extra': ['localIncidenceAngle', 'incidenceAngleFromEllipsoid',
-                             'scatteringArea', 'layoverShadowMask', 'gammaSigmaRatio'],
             'dem_resampling_method': 'BILINEAR_INTERPOLATION',
             'img_resampling_method': 'BILINEAR_INTERPOLATION',
             'slc_clean_edges': True,

diff --git a/S1_NRB/metadata/extract.py b/S1_NRB/metadata/extract.py
@@ -391,7 +391,8 @@ def _get_block_offset(band):
 def calc_geolocation_accuracy(swath_identifier, ei_tif, dem_type, etad):
     """
     Calculates the radial root mean square error, which is a target requirement of the CARD4L NRB specification
-    (Item 4.3). For more information see: https://s1-nrb.readthedocs.io/en/latest/general/geoaccuracy.html
+    (Item 4.3). For more information see: https://s1-nrb.readthedocs.io/en/latest/general/geoaccuracy.html.
+    Currently only the Copernicus DEM is supported.
     
     Parameters
     ----------
@@ -406,8 +407,8 @@ def calc_geolocation_accuracy(swath_identifier, ei_tif, dem_type, etad):
     
     Returns
     -------
-    rmse_planar: float
-        The calculated rRMSE value rounded to two decimal places.
+    rmse_planar: float or None
+        The calculated rRMSE value rounded to two decimal places or None if a DEM other than Copernicus is used.
     """
     if 'copernicus' not in dem_type.lower():
         return None
@@ -489,9 +490,13 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     sid0 = src_sid[list(src_sid.keys())[0]]  # first key/first file; used to extract some common metadata
     swath_id = re.search('_(IW|EW|S[1-6])_', os.path.basename(sid0.file)).group().replace('_', '')
 
-    ref_tif = finder(target, ['[hv]{2}-g-lin.tif$'], regex=True)[0]
+    ref_tif = finder(target, ['[hv]{2}-[gs]-lin.tif$'], regex=True)[0]
     np_tifs = finder(target, ['-np-[hv]{2}.tif$'], regex=True)
-    ei_tif = finder(target, ['-ei.tif$'], regex=True)[0]
+    ei_tifs = finder(target, ['-ei.tif$'], regex=True)
+    if len(ei_tifs) > 0:
+        ei_tif = ei_tifs[0]
+    else:
+        ei_tif = None
 
     product_id = os.path.basename(target)
     prod_meta = get_prod_meta(product_id=product_id, tif=ref_tif,
@@ -507,8 +512,11 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     tups = [(key, ITEM_MAP[key]['z_error']) for key in ITEM_MAP.keys()]
     z_err_dict = dict(tups)
 
-    geocorr_acc = calc_geolocation_accuracy(swath_identifier=swath_id, ei_tif=ei_tif,
-                                            dem_type=dem_type, etad=config['etad'])
+    if ei_tif is not None:
+        geocorr_acc = calc_geolocation_accuracy(swath_identifier=swath_id, ei_tif=ei_tif,
+                                                dem_type=dem_type, etad=config['etad'])
+    else:
+        geocorr_acc = None
 
     # Common metadata (sorted alphabetically)
     meta['common']['antennaLookDirection'] = 'RIGHT'
@@ -545,7 +553,7 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     meta['prod']['azimuthNumberOfLooks'] = prod_meta['ML_nAzLooks']
     meta['prod']['backscatterConvention'] = 'linear power'
     meta['prod']['backscatterConversionEq'] = '10*log10(DN)'
-    meta['prod']['backscatterMeasurement'] = 'gamma0'
+    meta['prod']['backscatterMeasurement'] = 'gamma0' if re.search('g-lin', ref_tif) else 'sigma0'
     meta['prod']['card4l-link'] = 'https://ceos.org/ard/files/PFS/NRB/v5.5/CARD4L-PFS_NRB_v5.5.pdf'
     meta['prod']['card4l-version'] = '5.5'
     meta['prod']['crsEPSG'] = str(prod_meta['epsg'])
@@ -586,7 +594,7 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     meta['prod']['griddingConvention'] = 'Military Grid Reference System (MGRS)'
     meta['prod']['licence'] = config['meta']['licence']
     meta['prod']['mgrsID'] = prod_meta['mgrsID']
-    meta['prod']['NRApplied'] = True if len(np_tifs) > 0 else False
+    meta['prod']['NRApplied'] = True
     meta['prod']['NRAlgorithm'] = 'https://sentinel.esa.int/documents/247904/2142675/Thermal-Denoising-of-Products-' \
                                   'Generated-by-Sentinel-1-IPF' if meta['prod']['NRApplied'] else None
     meta['prod']['numberOfAcquisitions'] = str(len(src_sid))

diff --git a/S1_NRB/metadata/mapping.py b/S1_NRB/metadata/mapping.py
@@ -17,6 +17,10 @@
             'vh-g-lin': {'z_error': 1e-4},
             'hh-g-lin': {'z_error': 1e-4},
             'hv-g-lin': {'z_error': 1e-4},
+            'vv-s-lin': {'z_error': 1e-4},
+            'vh-s-lin': {'z_error': 1e-4},
+            'hh-s-lin': {'z_error': 1e-4},
+            'hv-s-lin': {'z_error': 1e-4},
             'ei': {'z_error': 1e-3},
             'em': {'z_error': 1e-3},
             'dm': {'z_error': 0.0},
@@ -27,7 +31,8 @@
             'np-vv': {'z_error': 2e-5},
             'np-vh': {'z_error': 2e-5},
             'np-hh': {'z_error': 2e-5},
-            'np-hv': {'z_error': 2e-5}}
+            'np-hv': {'z_error': 2e-5},
+            'sg': {'z_error': 1e-4}}
 
 # Source data resolution
 # https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-1-sar/products-algorithms/level-1-algorithms/single-look-complex
@@ -140,7 +145,11 @@
               '-np-[vh]{2}.tif': {'type': 'Sigma-0',
                                   'unit': 'dB',
                                   'role': 'noise-power',
-                                  'title': 'Noise Power'}}
+                                  'title': 'Noise Power'},
+              '-sg.tif': {'type': 'Ratio',
+                          'unit': None,
+                          'role': 'sigma-gamma-ratio',
+                          'title': 'Sigma0 ellipsoidal to gamma0 RTC ratio'}}
 
 # https://sentinel.esa.int/documents/247904/1653442/Guide-to-Sentinel-1-Geocoding.pdf
 SLC_ACC_MAP = {'SM': {'ALE': {'rg': -3.02,

diff --git a/S1_NRB/metadata/stac.py b/S1_NRB/metadata/stac.py
@@ -190,15 +190,28 @@ def product_json(meta, target, tifs, exist_ok=False):
         size = os.path.getsize(asset)
 
         if 'measurement' in asset:
-            key = re.search('[a-z]{2}-[g|s]-(lin|log)', asset).group()
+            pattern = '(?P<key>(?P<pol>[vhc]{2})-(?P<nought>[gs])-(?P<scaling>lin|log))'
+            info = re.search(pattern, asset).groupdict()
+            key = info['key']
 
-            if key == 'cc-g-lin':
-                title = 'RGB color composite (vv-g-lin, vh-g-lin, vv-g-lin/vh-g-lin)'
+            if re.search('cc-[gs]-lin', key):
+                pols = meta['common']['polarisationChannels']
+                co = pols.pop(0) if pols[0][0] == pols[0][1] else pols.pop(1)
+                cross = pols[0]
+                title = 'RGB color composite (' \
+                        '{co}-{nought}-lin, ' \
+                        '{cross}-{nought}-lin, ' \
+                        '{co}-{nought}-lin/{cross}-{nought}-lin)'
+                title = title.format(co=co.lower(),
+                                     cross=cross.lower(),
+                                     nought=info['nought'])
             else:
-                pol = re.search('[vh]{2}', asset).group()
-                nought = measurement_title_dict[re.search('-[g|s]-', asset).group().replace('-', '')]
-                scaling = measurement_title_dict[re.search('(lin|log)', asset).group()]
-                title = '{} {} RTC backscatter, {} scaling'.format(pol.upper(), nought, scaling)
+                skeleton = '{pol} {nought} {subtype} backscatter, {scale} scaling'
+                subtype = 'RTC' if info['nought'] == 'g' else 'ellipsoidal'
+                title = skeleton.format(pol=info['pol'].upper(),
+                                        nought=measurement_title_dict[info['nought']],
+                                        subtype=subtype,
+                                        scale=measurement_title_dict[info['scaling']])
 
             header_size = get_header_size(tif=asset)
             if asset.endswith('.tif'):