diff --git a/S1_NRB/ancillary.py b/S1_NRB/ancillary.py index 87230763..9e238f83 100644 --- a/S1_NRB/ancillary.py +++ b/S1_NRB/ancillary.py @@ -193,7 +193,7 @@ def set_logging(config, debug=False): ---------- config: dict Dictionary of the parsed config parameters for the current process. - debug: bool, optional + debug: bool Set pyroSAR logging level to DEBUG? Default is False. Returns @@ -275,7 +275,7 @@ def _log_process_config(logger, config): Parameters ---------- - logger: Logger + logger: logging.Logger The logger to which the header is added to. config: dict Dictionary of the parsed config parameters for the current process. @@ -344,7 +344,7 @@ def log(handler, mode, proc_step, scenes, msg): scenes: str or list[str] Scenes that are currently being processed. msg: Any - The massage that should be logged. + The message that should be logged. """ proc_step = proc_step.zfill(7).replace('0', ' ') message = '[{proc_step}] -- {scenes} -- {msg}' diff --git a/S1_NRB/dem.py b/S1_NRB/dem.py index 62719e1b..76dbb952 100644 --- a/S1_NRB/dem.py +++ b/S1_NRB/dem.py @@ -221,11 +221,11 @@ def mosaic(geometry, dem_type, outname, epsg=None, kml_file=None, The DEM type. outname: str The name of the mosaic. - epsg: int + epsg: int or None The coordinate reference system as an EPSG code. - kml_file: str + kml_file: str or None The KML file containing the MGRS tile geometries. - dem_dir: str + dem_dir: str or None The directory containing the DEM MGRS tiles. username: str or None The username for accessing the DEM tiles. If None and authentication is required diff --git a/S1_NRB/metadata/extract.py b/S1_NRB/metadata/extract.py index e476d59c..5902873e 100644 --- a/S1_NRB/metadata/extract.py +++ b/S1_NRB/metadata/extract.py @@ -81,7 +81,7 @@ def vec_from_srccoords(coord_list): Parameters ---------- - coord_list: list[list[tuple(float, float)]] + coord_list: list[list[tuple[float]]] List containing for each source scene a list of coordinate pairs as retrieved from the metadata stored in an :class:`~pyroSAR.drivers.ID` object. @@ -126,13 +126,13 @@ def vec_from_srccoords(coord_list): def etree_from_sid(sid): """ - Retrieve the manifest and annotation XML data of a scene as a dictionary using an :class:`~pyroSAR.drivers.ID` + Retrieve the manifest and annotation XML data of a scene as a dictionary using an :class:`pyroSAR.drivers.ID` object. Parameters ---------- - sid: :class:`pyroSAR.drivers.ID` - A pyroSAR :class:`~pyroSAR.drivers.ID` object generated with :func:`pyroSAR.drivers.identify`. + sid: pyroSAR.drivers.ID + A pyroSAR :class:`~pyroSAR.drivers.ID` object generated with e.g. :func:`pyroSAR.drivers.identify`. Returns ------- @@ -163,13 +163,13 @@ def geometry_from_vec(vectorobject): Parameters ---------- - vectorobject: :class:`spatialist.vector.Vector` + vectorobject: spatialist.vector.Vector The vector object to extract geometry information from. Returns ------- out: dict - A dictionary containing the geometry information extracted from the vectorobject. + A dictionary containing the geometry information extracted from the vector object. """ out = {} vec = vectorobject @@ -206,21 +206,21 @@ def find_in_annotation(annotation_dict, pattern, single=False, out_type='str'): A dict of annotation files in the form: {'swath ID': `lxml.etree._Element` object} pattern: str The pattern to search for in each annotation file. - single: bool, optional + single: bool If True, the results found in each annotation file are expected to be the same and therefore only a single value will be returned instead of a dict. If the results differ, an error is raised. Default is False. - out_type: str, optional + out_type: str Output type to convert the results to. Can be one of the following: - - str (default) - - float - - int + - 'str' (default) + - 'float' + - 'int' Returns ------- out: dict A dictionary of the results containing a list for each of the annotation files. E.g., - {'swath ID': list[str, float or int]} + {'swath ID': list[str or float or int]} """ out = {} for s, a in annotation_dict.items(): @@ -305,10 +305,11 @@ def extract_pslr_islr(annotation_dict): Returns ------- - pslr: float - Mean PSLR value for all swaths of the scene. - islr: float - Mean ISLR value for all swaths of the scene. + tuple[float] + a tuple with the following values: + + - pslr: Mean PSLR value for all swaths of the scene. + - islr: Mean ISLR value for all swaths of the scene. """ swaths = list(annotation_dict.keys()) pslr_dict = find_in_annotation(annotation_dict=annotation_dict, pattern='.//crossCorrelationPslr', out_type='float') diff --git a/S1_NRB/metadata/stac.py b/S1_NRB/metadata/stac.py index 759d6ff2..e1f5fd40 100644 --- a/S1_NRB/metadata/stac.py +++ b/S1_NRB/metadata/stac.py @@ -27,7 +27,7 @@ def product_json(meta, target, tifs, exist_ok=False): A path pointing to the root directory of a product scene. tifs: list[str] List of paths to all GeoTIFF files of the currently processed NRB product. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ scene_id = os.path.basename(target) @@ -308,7 +308,7 @@ def source_json(meta, target, exist_ok=False): Metadata dictionary generated with :func:`~S1_NRB.metadata.extract.meta_dict`. target: str A path pointing to the root directory of a product scene. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ metadir = os.path.join(target, 'source') @@ -458,7 +458,7 @@ def parse(meta, target, tifs, exist_ok=False): A path pointing to the root directory of a product scene. tifs: list[str] List of paths to all GeoTIFF files of the currently processed NRB product. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ source_json(meta=meta, target=target, exist_ok=exist_ok) @@ -478,9 +478,9 @@ def make_catalog(directory, recursive=True, silent=False): ---------- directory: str Path to a directory that contains Sentinel-1 NRB products. - recursive: bool, optional + recursive: bool Search for NRB products in `directory` recursively? Default is True. - silent: bool, optional + silent: bool Should the output during directory reorganization be suppressed? Default is False. Returns @@ -574,11 +574,11 @@ def _reorganize_by_tile(directory, products=None, recursive=True, silent=False): ---------- directory: str Path to a directory that contains Sentinel-1 NRB products. - products: list[str], optional + products: list[str] or None List of NRB product paths. Will be created from `directory` if not provided. - recursive: bool, optional + recursive: bool Search for NRB products in `directory` recursively? Default is True. - silent: bool, optional + silent: bool If False (default), a message for each NRB product is printed if it has been moved to a new location or not. Returns diff --git a/S1_NRB/metadata/xml.py b/S1_NRB/metadata/xml.py index 012bba07..1ebf249a 100644 --- a/S1_NRB/metadata/xml.py +++ b/S1_NRB/metadata/xml.py @@ -20,7 +20,7 @@ def _om_time(root, nsmap, scene_id, time_start, time_stop): Parameters ---------- - root: etree.Element + root: lxml.etree.Element Root XML element. nsmap: dict Dictionary listing abbreviation (key) and URI (value) of all necessary XML namespaces. @@ -53,7 +53,7 @@ def _om_procedure(root, nsmap, scene_id, meta, uid=None, prod=True): Parameters ---------- - root: etree.Element + root: lxml.etree.Element Root XML element. nsmap: dict Dictionary listing abbreviation (key) and URI (value) of all necessary XML namespaces. @@ -61,9 +61,9 @@ def _om_procedure(root, nsmap, scene_id, meta, uid=None, prod=True): Scene basename. meta: dict Metadata dictionary generated with :func:`~S1_NRB.metadata.extract.meta_dict` - uid: str, optional + uid: str or None Unique identifier of a source SLC scene. - prod: bool, optional + prod: bool Return XML subelements for further usage in :func:`~S1_NRB.metadata.xml.product_xml` parsing function? Default is True. If False, the XML subelements for further usage in the :func:`~S1_NRB.metadata.xml.source_xml` parsing function will be returned. @@ -167,7 +167,7 @@ def _om_feature_of_interest(root, nsmap, scene_id, extent, center): Parameters ---------- - root: etree.Element + root: lxml.etree.Element Root XML element. nsmap: dict Dictionary listing abbreviation (key) and URI (value) of all necessary XML namespaces. @@ -213,7 +213,7 @@ def product_xml(meta, target, tifs, nsmap, exist_ok=False): List of paths to all GeoTIFF files of the currently processed NRB product. nsmap: dict Dictionary listing abbreviation (key) and URI (value) of all necessary XML namespaces. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ scene_id = os.path.basename(target) @@ -484,7 +484,7 @@ def source_xml(meta, target, nsmap, exist_ok=False): A path pointing to the root directory of a product scene. nsmap: dict Dictionary listing abbreviation (key) and URI (value) of all necessary XML namespaces. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ metadir = os.path.join(target, 'source') @@ -653,7 +653,7 @@ def parse(meta, target, tifs, exist_ok=False): A path pointing to the root directory of a product scene. tifs: list[str] List of paths to all GeoTIFF files of the currently processed NRB product. - exist_ok: bool, optional + exist_ok: bool Do not create files if they already exist? """ NS_MAP_prod = deepcopy(NS_MAP) diff --git a/S1_NRB/nrb.py b/S1_NRB/nrb.py index 4f391d83..8c9fda7f 100644 --- a/S1_NRB/nrb.py +++ b/S1_NRB/nrb.py @@ -47,15 +47,15 @@ def format(config, scenes, datadir, outdir, tile, extent, epsg, wbm=None, Spatial extent of the MGRS tile, derived from a :class:`~spatialist.vector.Vector` object. epsg: int The CRS used for the NRB product; provided as an EPSG code. - wbm: str, optional + wbm: str or None Path to a water body mask file with the dimensions of an MGRS tile. - multithread: bool, optional + 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`. - compress: str, optional + compress: str or None Compression algorithm to use. See https://gdal.org/drivers/raster/gtiff.html#creation-options for options. Defaults to 'LERC_DEFLATE'. - overviews: list[int], optional + overviews: list[int] or None Internal overview levels to be created for each GeoTIFF file. Defaults to [2, 4, 9, 18, 36] Returns @@ -350,22 +350,22 @@ def create_vrt(src, dst, fun, relpaths=False, scale=None, offset=None, args=None https://gdal.org/drivers/raster/vrt.html#default-pixel-functions. Furthermore, the option 'decibel' can be specified, which will implement a custom pixel function that uses Python code for decibel conversion (10*log10). - relpaths: bool, optional + relpaths: bool Should all `SourceFilename` XML elements with attribute `@relativeToVRT="0"` be updated to be paths relative to the output VRT file? Default is False. - scale: int, optional + scale: int or None The scale that should be applied when computing “real” pixel values from scaled pixel values on a raster band. Will be ignored if `fun='decibel'`. - offset: float, optional + offset: float or None The offset that should be applied when computing “real” pixel values from scaled pixel values on a raster band. Will be ignored if `fun='decibel'`. - args: dict, optional + args: dict or None arguments for `fun` passed as `PixelFunctionArguments`. Requires GDAL>=3.5 to be read. - options: dict, optional + options: dict or None Additional parameters passed to `gdal.BuildVRT`. - overviews: list[int], optional + overviews: list[int] or None Internal overview levels to be created for each raster file. - overview_resampling: str, optional + overview_resampling: str or None Resampling method for overview levels. Examples @@ -556,10 +556,10 @@ def calc_product_start_stop(src_ids, extent, epsg): Returns ------- - start: str - Start time of the NRB product formatted as `%Y%m%dT%H%M%S` in UTC. - stop: str - Stop time of the NRB product formatted as `%Y%m%dT%H%M%S` in UTC. + tuple[str] + + - Start time of the NRB product formatted as `%Y%m%dT%H%M%S` in UTC. + - Stop time of the NRB product formatted as `%Y%m%dT%H%M%S` in UTC. """ with bbox(extent, epsg) as tile_geom: tile_geom.reproject(4326) @@ -658,7 +658,7 @@ def create_data_mask(outname, datasets, extent, epsg, driver, creation_opt, Resampling method for overview levels. dst_nodata: int or str Nodata value to write to the output raster. - wbm: str, optional + wbm: str or None Path to a water body mask file with the dimensions of an MGRS tile. """ print(outname) diff --git a/S1_NRB/processor.py b/S1_NRB/processor.py index 03c82158..f665e262 100644 --- a/S1_NRB/processor.py +++ b/S1_NRB/processor.py @@ -8,6 +8,7 @@ from S1_NRB.config import get_config, snap_conf, gdal_conf import S1_NRB.ancillary as anc import S1_NRB.tile_extraction as tile_ex +from datetime import datetime, timedelta gdal.UseExceptions() @@ -90,7 +91,7 @@ def main(config_file, section_name='PROCESSING', debug=False): aoi_tiles = tile_ex.tile_from_aoi(vector=vec, kml=config['kml_file']) del vec #################################################################################################################### - # DEM download and WBM MGRS-tiling + # main SAR processing if rtc_flag: username, password = dem.authenticate(dem_type=config['dem_type'], username=None, password=None) for i, scene in enumerate(scenes): @@ -98,6 +99,17 @@ def main(config_file, section_name='PROCESSING', debug=False): out_dir_scene = os.path.join(config['rtc_dir'], scene_base) tmp_dir_scene = os.path.join(config['tmp_dir'], scene_base) + print(f'###### [ RTC] Scene {i + 1}/{len(scenes)}: {scene.scene}') + if os.path.isdir(out_dir_scene): + msg = 'Already processed - Skip!' + print('### ' + msg) + anc.log(handler=logger, mode='info', proc_step='GEOCODE', scenes=scene.scene, msg=msg) + continue + else: + os.makedirs(out_dir_scene) + os.makedirs(tmp_dir_scene, exist_ok=True) + ############################################################################################################ + # Preparation of DEM and WBM dem.prepare(vector=scene.bbox(), threads=gdal_prms['threads'], dem_dir=None, wbm_dir=config['wbm_dir'], dem_type=config['dem_type'], kml_file=config['kml_file'], @@ -116,24 +128,14 @@ def main(config_file, section_name='PROCESSING', debug=False): dem.mosaic(geometry=geom, outname=fname_dem, dem_type=config['dem_type'], username=username, password=password) ############################################################################################################ - # RTC processing - print(f'###### [ RTC] Scene {i + 1}/{len(scenes)}: {scene.scene}') - if os.path.isdir(out_dir_scene): - msg = 'Already processed - Skip!' - print('### ' + msg) - anc.log(handler=logger, mode='info', proc_step='GEOCODE', scenes=scene.scene, msg=msg) - continue - else: - os.makedirs(out_dir_scene) - os.makedirs(tmp_dir_scene, exist_ok=True) - ############################################### # ETAD correction if config['etad']: msg = '###### [ ETAD] Scene {s}/{s_total}: {scene}' print(msg.format(s=i + 1, s_total=len(scenes), scene=scene.scene)) scene = etad.process(scene=scene, etad_dir=config['etad_dir'], out_dir=tmp_dir_scene, log=logger) - ############################################### + ############################################################################################################ + # determination of look factors if scene.product == 'SLC': rlks = {'IW': 5, 'SM': 6, @@ -145,12 +147,29 @@ def main(config_file, section_name='PROCESSING', debug=False): azlks *= int(geocode_prms['spacing'] / 10) else: rlks = azlks = None - ############################################### + ############################################################################################################ + # get neighbouring GRD scenes to add a buffer to the geocoded scenes + # otherwise there will be a gap between final geocoded images. + neighbors = None + if scene.product == 'GRD': + f = '%Y%m%dT%H%M%S' + td = timedelta(seconds=2) + start = datetime.strptime(scene.start, f) - td + start = datetime.strftime(start, f) + stop = datetime.strptime(scene.stop, f) + td + stop = datetime.strftime(stop, f) + with Archive(config['db_file']) as db: + neighbors = db.select(mindate=start, maxdate=stop, date_strict=False, + sensor=scene.sensor, product=scene.product, + acquisition_mode=scene.acquisition_mode) + del neighbors[neighbors.index(scene.scene)] + ############################################################################################################ + # main processing routine start_time = time.time() try: snap.process(scene=scene.scene, outdir=config['rtc_dir'], tmpdir=config['tmp_dir'], kml=config['kml_file'], - dem=fname_dem, + dem=fname_dem, neighbors=neighbors, 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) @@ -173,6 +192,8 @@ def main(config_file, section_name='PROCESSING', debug=False): del vec # filter the tile selection based on the user geometry config tiles = [x for x in tiles if x.mgrs in aoi_tiles] + t_total = len(tiles) + s_total = len(selection_grouped) for t, tile in enumerate(tiles): outdir = os.path.join(config['nrb_dir'], tile.mgrs) os.makedirs(outdir, exist_ok=True) @@ -181,10 +202,10 @@ def main(config_file, section_name='PROCESSING', debug=False): wbm = None extent = tile.extent epsg = tile.getProjection('epsg') - msg = '###### [ NRB] Tile {t}/{t_total}: {tile} | Scenes {s}/{s_total}: {scenes} ' - print(msg.format(tile=tile.mgrs, t=t + 1, t_total=len(aoi_tiles), + 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], - s=s + 1, s_total=len(selection_grouped))) + s=s + 1, s_total=s_total)) try: msg = nrb.format(config=config, scenes=scenes_fnames, datadir=config['rtc_dir'], outdir=outdir, tile=tile.mgrs, extent=extent, epsg=epsg, diff --git a/S1_NRB/snap.py b/S1_NRB/snap.py index 827dd457..8b9d56c1 100644 --- a/S1_NRB/snap.py +++ b/S1_NRB/snap.py @@ -13,10 +13,9 @@ from S1_NRB.ancillary import get_max_ext -def mli(src, dst, workflow, spacing=None, rlks=None, azlks=None, allow_res_osv=True): +def mli(src, dst, workflow, spacing=None, rlks=None, azlks=None): """ - Create a multi-looked image (MLI). The following operators are used (optional steps in brackets): - Apply-Orbit-File(->Remove-GRD-Border-Noise)->Calibration->ThermalNoiseRemoval(->TOPSAR-Deburst->Multilook) + Multi-looing. Parameters ---------- @@ -30,10 +29,48 @@ def mli(src, dst, workflow, spacing=None, rlks=None, azlks=None, allow_res_osv=T the target pixel spacing for automatic determination of looks using function :func:`pyroSAR.ancillary.multilook_factors`. Overridden by arguments `rlks` and `azlks` if they are not None. - rlks: int + rlks: int or None the number of range looks. - azlks: int + azlks: int or None the number of azimuth looks. + + Returns + ------- + + """ + scene = identify(src) + wf = parse_recipe('blank') + ############################################ + read = parse_node('Read') + read.parameters['file'] = scene.scene + wf.insert_node(read) + ############################################ + ml = mli_parametrize(scene=scene, spacing=spacing, rlks=rlks, azlks=azlks) + if ml is not None: + wf.insert_node(ml, before=read.id) + ############################################ + write = parse_node('Write') + wf.insert_node(write, before=ml.id) + write.parameters['file'] = dst + write.parameters['formatName'] = 'BEAM-DIMAP' + ############################################ + wf.write(workflow) + gpt(xmlfile=workflow, tmpdir=os.path.dirname(dst)) + + +def pre(src, dst, workflow, allow_res_osv=True): + """ + SAR preprocessing. The following operators are used (optional steps in brackets): + Apply-Orbit-File(->Remove-GRD-Border-Noise)->Calibration->ThermalNoiseRemoval(->TOPSAR-Deburst) + + Parameters + ---------- + src: str + the file name of the source scene + dst: str + the file name of the target scene. Format is BEAM-DIMAP. + workflow: str + the output SNAP XML workflow filename. allow_res_osv: bool Also allow the less accurate RES orbit files to be used? @@ -77,11 +114,6 @@ def mli(src, dst, workflow, spacing=None, rlks=None, azlks=None, allow_res_osv=T wf.insert_node(deb, before=last.id) last = deb ############################################ - ml = mli_parametrize(scene=scene, spacing=spacing, rlks=rlks, azlks=azlks) - if ml is not None: - wf.insert_node(ml, before=last.id) - last = ml - ############################################ write = parse_node('Write') wf.insert_node(write, before=last.id) write.parameters['file'] = dst @@ -91,7 +123,65 @@ def mli(src, dst, workflow, spacing=None, rlks=None, azlks=None, allow_res_osv=T gpt(xmlfile=workflow, tmpdir=os.path.dirname(dst)) -def rtc(src, dst, workflow, dem, dem_resampling_method='BILINEAR_INTERPOLATION', sigma0=True, scattering_area=True): +def grd_buffer(src, dst, workflow, neighbors, buffer=10): + """ + GRDs, unlike SLCs, do not overlap in azimuth. + With this function, a GRD can be buffered using the neighboring acquisitions. + First, all images are mosaicked using the `SliceAssembly` operator + and then subsetted to the extent of the main scene including a buffer. + + Parameters + ---------- + src: str + the file name of the source scene + dst: str + the file name of the target scene. Format is BEAM-DIMAP. + workflow: str + the output SNAP XML workflow filename. + neighbors: list[str] + the file names of neighboring scenes + buffer: int + the number of pixels to buffer + + Returns + ------- + + """ + scenes = identify_many([src] + neighbors, sortkey='start') + wf = parse_recipe('blank') + ############################################ + read_ids = [] + for scene in scenes: + read = parse_node('Read') + read.parameters['file'] = scene.scene + wf.insert_node(read) + read_ids.append(read.id) + ############################################ + asm = parse_node('SliceAssembly') + wf.insert_node(asm, before=read_ids) + ############################################ + id_main = [x.scene for x in scenes].index(src) + xmin = 0 + width = scenes[id_main].samples + ymin = 0 if id_main == 0 else scenes[0].lines - buffer + height = scenes[id_main].lines + buffer * 2 + sub = parse_node('Subset') + sub.parameters['region'] = [xmin, ymin, width, height] + sub.parameters['geoRegion'] = '' + sub.parameters['copyMetadata'] = True + wf.insert_node(sub, before=asm.id) + ############################################ + write = parse_node('Write') + wf.insert_node(write, before=sub.id) + write.parameters['file'] = dst + write.parameters['formatName'] = 'BEAM-DIMAP' + ############################################ + wf.write(workflow) + gpt(xmlfile=workflow, tmpdir=os.path.dirname(dst)) + + +def rtc(src, dst, workflow, dem, dem_resampling_method='BILINEAR_INTERPOLATION', + sigma0=True, scattering_area=True): """ Radiometric Terrain Flattening. @@ -109,7 +199,7 @@ def rtc(src, dst, workflow, dem, dem_resampling_method='BILINEAR_INTERPOLATION', the DEM resampling method. sigma0: bool output sigma0 RTC backscatter? - scattering_area:bool + scattering_area: bool output scattering area image? Returns @@ -218,7 +308,7 @@ def geo(*src, dst, workflow, spacing, crs, geometry=None, buffer=0.01, Parameters ---------- - src + src: list[str or None] variable number of input scene file names dst: str the file name of the target scene. Format is BEAM-DIMAP. @@ -251,7 +341,7 @@ def geo(*src, dst, workflow, spacing, crs, geometry=None, buffer=0.01, img_resampling_method: str the SAR image resampling method bands - band ids for the input scenes in `args` as lists with keys bands, + band ids for the input scenes in `src` as lists with keys bands, e.g., bands1=['NESZ_VV'], bands2=['Gamma0_VV'], ... Returns @@ -309,7 +399,7 @@ def process(scene, outdir, spacing, kml, dem, img_resampling_method='BILINEAR_INTERPOLATION', rlks=None, azlks=None, tmpdir=None, export_extra=None, allow_res_osv=True, slc_clean_edges=True, slc_clean_edges_pixels=4, - cleanup=True): + neighbors=None, cleanup=True): """ Main function for RTC processing with SNAP. @@ -352,6 +442,10 @@ def process(scene, outdir, spacing, kml, dem, Does not apply to layover-shadow mask. slc_clean_edges_pixels: int The number of pixels to erode. + neighbors: list[str] or None + (only applies to GRD) an optional list of neighboring scenes to add + a buffer around the main scene. If GRDs are processed compeletely + independently, gaps are introduced due to a missing overlap between GRDs. cleanup: bool Delete intermediate files after successful process termination? @@ -388,13 +482,49 @@ def process(scene, outdir, spacing, kml, dem, workflows = [] ############################################################################ # general pre-processing + out_pre = tmp_base + '_pre.dim' + out_pre_wf = out_pre.replace('.dim', '.xml') + workflows.append(out_pre_wf) + if not os.path.isfile(out_pre): + pre(src=scene, dst=out_pre, workflow=out_pre_wf, + allow_res_osv=allow_res_osv) + ############################################################################ + # GRD buffering + if neighbors is not None: + # general preprocessing of neighboring scenes + out_pre_neighbors = [] + for item in neighbors: + basename_nb = os.path.splitext(os.path.basename(item))[0] + tmpdir_nb = os.path.join(tmpdir, basename_nb) + os.makedirs(tmpdir_nb, exist_ok=True) + tmp_base_nb = os.path.join(tmpdir_nb, basename_nb) + out_pre_nb = tmp_base_nb + '_pre.dim' + out_pre_nb_wf = out_pre_nb.replace('.dim', '.xml') + if not os.path.isfile(out_pre_nb): + print('### preprocessing neighbor:', item) + pre(src=item, dst=out_pre_nb, workflow=out_pre_nb_wf, + allow_res_osv=allow_res_osv) + out_pre_neighbors.append(out_pre_nb) + ######################################################################## + # buffering + out_buffer = tmp_base + '_buf.dim' + out_buffer_wf = out_buffer.replace('.dim', '.xml') + workflows.append(out_buffer_wf) + if not os.path.isfile(out_buffer): + grd_buffer(src=out_pre, dst=out_buffer, workflow=out_buffer_wf, + neighbors=out_pre_neighbors) + out_pre = out_buffer + ############################################################################ + # multi-looking out_mli = tmp_base + '_mli.dim' out_mli_wf = out_mli.replace('.dim', '.xml') - workflows.append(out_mli_wf) if not os.path.isfile(out_mli): - mli(src=scene, dst=out_mli, workflow=out_mli_wf, - spacing=spacing, rlks=rlks, azlks=azlks, - allow_res_osv=allow_res_osv) + mli(src=out_pre, dst=out_mli, workflow=out_mli_wf, + spacing=spacing, rlks=rlks, azlks=azlks) + if not os.path.isfile(out_mli): + out_mli = out_pre + else: + workflows.append(out_mli_wf) ############################################################################ # radiometric terrain flattening out_rtc = tmp_base + '_rtc.dim' @@ -523,22 +653,22 @@ def find_datasets(scene, outdir, epsg): subdir = os.path.join(scenedir, basename + f'_geo_{epsg}.data') if not os.path.isdir(subdir): return - lookup = {'dm': r'layoverShadowMask\.img', - 'ei': r'incidenceAngleFromEllipsoid\.img', - 'gs': r'gammaSigmaRatio_[VH]{2}\.img', - 'lc': r'simulatedImage_[VH]{2}\.img', - 'li': r'localIncidenceAngle\.img'} + lookup = {'dm': r'layoverShadowMask\.img$', + 'ei': r'incidenceAngleFromEllipsoid\.img$', + 'gs': r'gammaSigmaRatio_[VH]{2}\.img$', + 'lc': r'simulatedImage_[VH]{2}\.img$', + 'li': r'localIncidenceAngle\.img$'} out = {} for key, pattern in lookup.items(): match = finder(target=subdir, matchlist=[pattern], regex=True) if len(match) > 0: out[key] = match[0] - pattern = r'Gamma0_(?P[VH]{2})\.img' + pattern = r'Gamma0_(?P[VH]{2})\.img$' gamma = finder(target=subdir, matchlist=[pattern], regex=True) for item in gamma: pol = re.search(pattern, item).group('pol') out[f'{pol.lower()}-g-lin'] = item - pattern = r'NESZ_(?P[VH]{2})\.img' + pattern = r'NESZ_(?P[VH]{2})\.img$' nesz = finder(target=subdir, matchlist=[pattern], regex=True) for item in nesz: pol = re.search(pattern, item).group('pol') @@ -564,12 +694,12 @@ def get_metadata(scene, outdir): """ basename = os.path.splitext(os.path.basename(scene))[0] scenedir = os.path.join(outdir, basename) - wf_mli = finder(scenedir, ['*mli.xml'])[0] - wf = parse_recipe(wf_mli) - if 'Multilook' in wf.operators: - rlks = int(wf['Multilook'].parameters['nRgLooks']) - azlks = int(wf['Multilook'].parameters['nAzLooks']) - else: - rlks = azlks = 1 + rlks = azlks = 1 + wf_mli = finder(scenedir, ['*mli.xml']) + if len(wf_mli) > 0: + wf = parse_recipe(wf_mli) + if 'Multilook' in wf.operators: + rlks = int(wf['Multilook'].parameters['nRgLooks']) + azlks = int(wf['Multilook'].parameters['nAzLooks']) return {'azlks': azlks, 'rlks': rlks} diff --git a/S1_NRB/tile_extraction.py b/S1_NRB/tile_extraction.py index 73ffa17b..6fac044d 100644 --- a/S1_NRB/tile_extraction.py +++ b/S1_NRB/tile_extraction.py @@ -48,6 +48,7 @@ def tile_from_aoi(vector, kml, epsg=None, strict=True, return_geometries=False): if return_geometries: sortkey = lambda x: x.mgrs with Vector(kml, driver='KML') as vec: + tilenames_src = [] tiles = [] for vector in vectors: vector.layer.ResetReading() @@ -56,7 +57,8 @@ def tile_from_aoi(vector, kml, epsg=None, strict=True, return_geometries=False): vec.layer.SetSpatialFilter(geom) for tile in vec.layer: tilename = tile.GetField('Name') - if tilename not in tiles: + if tilename not in tilenames_src: + tilenames_src.append(tilename) attrib = description2dict(tile.GetField('Description')) reproject = False if epsg is not None and attrib['EPSG'] not in epsg: