Merge pull request #85 from SAR-ARD/feature/snap_proc

general processor improvements

S1_NRB/ancillary.py

@@ -318,6 +318,7 @@ def _log_process_config(logger, config):
     db_file             {config['db_file']}
     kml_file            {config['kml_file']}
     gdal_threads        {config.get('gdal_threads')}
+    snap_gpt_args       {config['snap_gpt_args']}
     
     ====================================================================================================================
     SOFTWARE

S1_NRB/cli.py

@@ -26,6 +26,10 @@ def cli(ctx, config_file, section, debug, version):
     
     s1_nrb -c config.ini --acq_mode IW --annotation dm,id
     
+    The snap_gpt_args argument can be provided by using a dedicated argument separator and quotes:
+    
+    s1_nrb -c config.ini -- --snap_gpt_args "-J-Xmx100G -c 75G -q 30"
+    
     \b
     The following defaults are set:
     (processing section)
@@ -36,6 +40,7 @@ def cli(ctx, config_file, section, debug, version):
     - etad:              False
     - etad_dir:          None
     - gdal_threads:      4
+    - snap_gpt_args:     None
     - log_dir:           LOG
     - measurement:       gamma
     - nrb_dir:           NRB

S1_NRB/config.py

@@ -25,7 +25,7 @@ def get_keys(section):
                 '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', 'annotation', 'stac_catalog', 'stac_collections',
-                'sensor', 'date_strict']
+                'sensor', 'date_strict', 'snap_gpt_args']
     elif section == 'metadata':
         return ['access_url', 'licence', 'doi', 'processing_center']
     else:
@@ -51,7 +51,8 @@ def get_config(config_file, proc_section='PROCESSING', **kwargs):
                                        converters={'_annotation': _parse_annotation,
                                                    '_datetime': _parse_datetime,
                                                    '_stac_collections': _parse_list,
-                                                   '_tile_list': _parse_tile_list})
+                                                   '_tile_list': _parse_tile_list,
+                                                   '_list': _parse_list})
     if isinstance(config_file, str):
         if not os.path.isfile(config_file):
             raise FileNotFoundError("Config file {} does not exist.".format(config_file))
@@ -87,8 +88,10 @@ def get_config(config_file, proc_section='PROCESSING', **kwargs):
         proc_sec['dem_type'] = 'Copernicus 30m Global DEM'
     if 'date_strict' not in proc_sec.keys():
         proc_sec['date_strict'] = 'True'
-
-        # use previous defaults for measurement and annotation if they have not been defined
+    if 'snap_gpt_args' not in proc_sec.keys():
+        proc_sec['snap_gpt_args'] = 'None'
+
+    # use previous defaults for measurement and annotation if they have not been defined
     if 'measurement' not in proc_sec.keys():
         proc_sec['measurement'] = 'gamma'
     if 'annotation' not in proc_sec.keys():
@@ -173,6 +176,8 @@ def get_config(config_file, proc_section='PROCESSING', **kwargs):
             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)
+        if k == 'snap_gpt_args':
+            v = proc_sec.get_list(k)
         out_dict[k] = v
 
     if out_dict['db_file'] is None and out_dict['stac_catalog'] is None:
@@ -243,7 +248,10 @@ def _parse_list(s):
     if s in ['', 'None']:
         return None
     else:
-        return s.replace(' ', '').split(',')
+        if re.search(',', s):
+            return s.replace(' ', '').split(',')
+        else:
+            return s.split()
 
 
 def _keyval_check(key, val, allowed_keys):
@@ -279,8 +287,8 @@ def snap_conf(config):
             'allow_res_osv': True,
             'dem_resampling_method': 'BILINEAR_INTERPOLATION',
             'img_resampling_method': 'BILINEAR_INTERPOLATION',
-            'slc_clean_edges': True,
-            'slc_clean_edges_pixels': 4,
+            'clean_edges': True,
+            'clean_edges_pixels': 4,
             'cleanup': True
             }
 

S1_NRB/dem.py

@@ -302,7 +302,7 @@ def mosaic(geometry, dem_type, outname, epsg=None, kml_file=None,
             dem_autoload([geometry], demType=dem_type,
                          vrt=vrt, buffer=buffer, product='dem',
                          username=username, password=password)
-            dem_create(src=vrt, dst=outname, pbar=True,
+            dem_create(src=vrt, dst=outname, pbar=False,
                        geoid_convert=geoid_convert, geoid=geoid,
                        threads=threads, nodata=-32767)
 

S1_NRB/metadata/extract.py

@@ -279,7 +279,7 @@ def calc_performance_estimates(files):
     """
     out = {}
     for f in files:
-        pol = re.search('[vh]{2}', f).group().upper()
+        pol = re.search('np-([vh]{2})', f).group(1).upper()
         with Raster(f) as ras:
             arr = ras.array()
             # The following need to be of type float, not numpy.float32 in order to be JSON serializable
@@ -428,6 +428,10 @@ def calc_geolocation_accuracy(swath_identifier, ei_tif, dem_type, etad):
         stats = ras.allstats(approximate=False)
         ei_min = stats[0]['min']
 
+    if ei_min == 0:
+        raise RuntimeError(f'minimum ellipsoid incidence angle cannot be 0\n'
+                           f'(file: {ei_tif})')
+
     # Remove generated '.aux.xml' file
     aux = finder(os.path.dirname(ei_tif), ['.tif.aux.xml$'], regex=True, recursive=False)
     for file in aux:

S1_NRB/nrb.py

@@ -39,7 +39,9 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None,
     config: dict
         Dictionary of the parsed config parameters for the current process.
     scenes: list[str]
-        List of scenes to process. Either an individual scene or multiple, matching scenes (consecutive acquisitions).
+        List of scenes to process. Either a single scene or multiple, matching scenes (consecutive acquisitions).
+        All scenes are expected to overlap with `extent` and an error will be thrown if the processing output
+        cannot be found for any of the scenes.
     datadir: str
         The directory containing the datasets processed from the source scenes using pyroSAR.
     outdir: str
@@ -54,7 +56,7 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None,
         Path to a water body mask file with the dimensions of an MGRS tile.
     dem_type: str or None
         if defined, a DEM layer will be added to the product. The suffix `em` (elevation model) is used.
-        Default `None`: do not add a DEm layer.
+        Default `None`: do not add a DEM layer.
     multithread: bool
         Should `gdalwarp` use multithreading? Default is True. The number of threads used, can be adjusted in the
         `config.ini` file with the parameter `gdal_threads`.
@@ -105,8 +107,8 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None,
 
     src_ids, datasets = get_datasets(scenes=scenes, datadir=datadir, extent=extent, epsg=epsg)
     if len(src_ids) == 0:
-        raise RuntimeError('None of the scenes overlap with the current tile {tile_id}: '
-                           '\n{scenes}'.format(tile_id=tile, scenes=scenes))
+        print('None of the processed scenes overlap with the current tile {tile_id}'.format(tile_id=tile))
+        return
 
     if annotation is not None:
         allowed = []
@@ -345,16 +347,25 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None,
 
 def get_datasets(scenes, datadir, extent, epsg):
     """
+    Collect processing output for a list of scenes.
     Reads metadata from all source SLC/GRD scenes, finds matching output files in `datadir`
-    and filters both lists depending on the actual overlap of each SLC/GRD footprint
-    with the current MGRS tile geometry.
+    and filters both lists depending on the actual overlap of each SLC/GRD valid data coverage
+    with the current MGRS tile geometry. If no output is found for any scene the function will raise an error.
+    To obtain the extent of valid data coverage, first a binary
+    mask raster file is created with the name `datamask.tif`, which is stored in the same folder as
+    the processing output as found by :func:`~S1_NRB.snap.find_datasets`. Then, the boundary of this
+    binary mask is computed and stored as `datamask.gpkg` (see function :func:`spatialist.vector.boundary`).
+    If the provided `extent` does not overlap with this boundary, the output is discarded. This scenario
+    might occur when the scene's geometry read from its metadata overlaps with the tile but the actual
+    extent of data does not.
 
     Parameters
     ----------
     scenes: list[str]
         List of scenes to process. Either an individual scene or multiple, matching scenes (consecutive acquisitions).
     datadir: str
         The directory containing the datasets processed from the source scenes using pyroSAR.
+        The function will raise an error if the processing output cannot be found for all scenes in `datadir`.
     extent: dict
         Spatial extent of the MGRS tile, derived from a :class:`~spatialist.vector.Vector` object.
     epsg: int
@@ -367,15 +378,20 @@ def get_datasets(scenes, datadir, extent, epsg):
     datasets: list[dict]
         List of RTC processing output files that match each :class:`~pyroSAR.drivers.ID` object of `ids`.
         The format is a list of dictionaries per scene with keys as described by e.g. :func:`S1_NRB.snap.find_datasets`.
+    
+    See Also
+    --------
+    :func:`S1_NRB.snap.find_datasets`
     """
     ids = identify_many(scenes)
     datasets = []
-    for _id in ids:
+    for i, _id in enumerate(ids):
         files = find_datasets(scene=_id.scene, outdir=datadir, epsg=epsg)
         if files is not None:
             datasets.append(files)
-    if len(datasets) == 0:
-        raise RuntimeError("No datasets were found in the directory '{}'".format(datadir))
+        else:
+            base = os.path.basename(_id.scene)
+            raise RuntimeError(f'cannot find processing output for scene {base}')
 
     i = 0
     while i < len(datasets):
@@ -387,19 +403,21 @@ def get_datasets(scenes, datadir, extent, epsg):
             with Raster(measurements[0]) as ras:
                 arr = ras.array()
                 mask = ~np.isnan(arr)
+                # remove scene if file does not contain valid data
+                if len(mask[mask == 1]) == 0:
+                    del ids[i]
+                    del datasets[i]
+                    continue
                 with vectorize(target=mask, reference=ras) as vec:
                     with boundary(vec, expression="value=1") as bounds:
                         if not os.path.isfile(dm_ras):
-                            print('creating raster mask', i)
                             rasterize(vectorobject=bounds, reference=ras, outname=dm_ras)
                         if not os.path.isfile(dm_vec):
-                            print('creating vector mask', i)
                             bounds.write(outfile=dm_vec)
         with Vector(dm_vec) as bounds:
             with bbox(extent, epsg) as tile_geom:
                 inter = intersect(bounds, tile_geom)
                 if inter is None:
-                    print('removing dataset', i)
                     del ids[i]
                     del datasets[i]
                 else:

S1_NRB/processor.py

@@ -1,7 +1,7 @@
 import os
 import time
 from osgeo import gdal
-from spatialist import Vector, bbox
+from spatialist import Vector, bbox, intersect
 from spatialist.ancillary import finder
 from pyroSAR import identify_many, Archive
 from S1_NRB import etad, dem, nrb, snap
@@ -55,7 +55,7 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
     else:
         acq_mode_search = config['acq_mode']
 
-    vec = [None]
+    vec = None
     aoi_tiles = None
     selection = []
     if config['aoi_tiles'] is not None:
@@ -76,28 +76,53 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
         archive = STACArchive(url=config['stac_catalog'],
                               collections=config['stac_collections'])
 
+    # derive geometries and tiles from scene footprints
+    if vec is None:
+        selection_tmp = archive.select(sensor=config['sensor'],
+                                       vectorobject=None,
+                                       product=config['product'],
+                                       acquisition_mode=acq_mode_search,
+                                       mindate=config['mindate'],
+                                       maxdate=config['maxdate'],
+                                       date_strict=config['date_strict'])
+        scenes = identify_many(scenes=selection_tmp)
+        scenes_geom = [x.geometry() for x in scenes]
+        # select all tiles overlapping with the scenes for further processing
+        vec = tile_ex.tile_from_aoi(vector=scenes_geom, kml=config['kml_file'],
+                                    return_geometries=True)
+        aoi_tiles = [x.mgrs for x in vec]
+        del scenes_geom, scenes
+        # extend the time range to fully cover all tiles
+        # (one additional scene needed before and after each data take group)
+        mindate = config['mindate'] - timedelta(minutes=1)
+        maxdate = config['maxdate'] + timedelta(minutes=1)
+    else:
+        mindate = config['mindate']
+        maxdate = config['maxdate']
+
     for item in vec:
         selection.extend(
             archive.select(sensor=config['sensor'],
                            vectorobject=item,
                            product=config['product'],
                            acquisition_mode=acq_mode_search,
-                           mindate=config['mindate'],
-                           maxdate=config['maxdate'],
+                           mindate=mindate,
+                           maxdate=maxdate,
                            date_strict=config['date_strict']))
     selection = list(set(selection))
     del vec
 
     if len(selection) == 0:
         message = "No scenes could be found for the following search query:\n" \
+                  " sensor:    '{sensor}'\n" \
                   " product:   '{product}'\n" \
                   " acq. mode: '{acq_mode}'\n" \
                   " mindate:   '{mindate}'\n" \
                   " maxdate:   '{maxdate}'\n"
         raise RuntimeError(message.format(acq_mode=config['acq_mode'], product=config['product'],
                                           mindate=config['mindate'], maxdate=config['maxdate'],
                                           scene_dir=config['scene_dir']))
-    scenes = identify_many(selection)
+    scenes = identify_many(selection, sortkey='start')
     anc.check_acquisition_completeness(scenes=scenes, archive=archive)
     archive.close()
 
@@ -206,6 +231,7 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
                              tmpdir=config['tmp_dir'], kml=config['kml_file'],
                              dem=fname_dem, neighbors=neighbors,
                              export_extra=export_extra,
+                             gpt_args=config['snap_gpt_args'],
                              rlks=rlks, azlks=azlks, **geocode_prms)
                 t = round((time.time() - start_time), 2)
                 anc.log(handler=logger, mode='info', proc_step='RTC', scenes=scene.scene, msg=t)
@@ -215,7 +241,7 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
     ####################################################################################################################
     # NRB - final product generation
     if nrb_flag:
-        # prepare DEM and WBM MGRS tiles
+        # prepare WBM MGRS tiles
         vec = [x.geometry() for x in scenes]
         extent = anc.get_max_ext(geometries=vec)
         del vec
@@ -225,10 +251,9 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
                         dem_type=config['dem_type'], kml_file=config['kml_file'],
                         tilenames=aoi_tiles, username=username, password=password,
                         dem_strict=True)
-
+        print('preparing NRB products')
         selection_grouped = anc.group_by_time(scenes=scenes)
         for s, scenes in enumerate(selection_grouped):
-            scenes_fnames = [x.scene for x in scenes]
             # check that the scenes can really be grouped together
             anc.check_scene_consistency(scenes=scenes)
             # get the geometries of all tiles that overlap with the current scene group
@@ -241,6 +266,9 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
             t_total = len(tiles)
             s_total = len(selection_grouped)
             for t, tile in enumerate(tiles):
+                # select all scenes from the group whose footprint overlaps with the current tile
+                scenes_sub = [x for x in scenes if intersect(tile, x.geometry())]
+                scenes_sub_fnames = [x.scene for x in scenes_sub]
                 outdir = os.path.join(config['nrb_dir'], tile.mgrs)
                 os.makedirs(outdir, exist_ok=True)
                 fname_wbm = os.path.join(config['wbm_dir'], config['dem_type'],
@@ -253,19 +281,19 @@ def main(config_file, section_name='PROCESSING', debug=False, **kwargs):
                 epsg = tile.getProjection('epsg')
                 msg = '###### [    NRB] Tile {t}/{t_total}: {tile} | Scenes: {scenes} '
                 print(msg.format(tile=tile.mgrs, t=t + 1, t_total=t_total,
-                                 scenes=[os.path.basename(s.scene) for s in scenes],
+                                 scenes=[os.path.basename(s) for s in scenes_sub_fnames],
                                  s=s + 1, s_total=s_total))
                 try:
-                    msg = nrb.format(config=config, scenes=scenes_fnames, datadir=config['rtc_dir'],
+                    msg = nrb.format(config=config, scenes=scenes_sub_fnames, datadir=config['rtc_dir'],
                                      outdir=outdir, tile=tile.mgrs, extent=extent, epsg=epsg,
                                      wbm=fname_wbm, dem_type=nrb_dem_type, kml=config['kml_file'],
                                      multithread=gdal_prms['multithread'], annotation=config['annotation'],
                                      update=update)
                     if msg == 'Already processed - Skip!':
                         print('### ' + msg)
-                    anc.log(handler=logger, mode='info', proc_step='NRB', scenes=scenes_fnames, msg=msg)
+                    anc.log(handler=logger, mode='info', proc_step='NRB', scenes=scenes_sub_fnames, msg=msg)
                 except Exception as e:
-                    anc.log(handler=logger, mode='exception', proc_step='NRB', scenes=scenes_fnames, msg=e)
+                    anc.log(handler=logger, mode='exception', proc_step='NRB', scenes=scenes_sub_fnames, msg=e)
                     continue
             del tiles
         gdal.SetConfigOption('GDAL_NUM_THREADS', gdal_prms['threads_before'])