From 24b6fb7ee5041bf4136b465c5138652d89bef079 Mon Sep 17 00:00:00 2001 From: Guido Riembauer Date: Tue, 16 Apr 2024 14:20:43 +0200 Subject: [PATCH 1/2] first_version_aoi_from_region --- .../i.sentinel1.pyrosargeocode.html | 6 +- .../i.sentinel1.pyrosargeocode.py | 63 +++++++++++++++++-- 2 files changed, 61 insertions(+), 8 deletions(-) diff --git a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.html b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.html index f5f45f44..3ddf0022 100644 --- a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.html +++ b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.html @@ -30,8 +30,10 @@

DESCRIPTION

the l-flag, m-flag, or r-flag are given.

-If only a specifc extent of the Sentinel-1 images is supposed to be geocoded, the user can define -the Area of Interest aoi in form of a GeoJSON with a single geometry. +If only a specific extent of the Sentinel-1 images is supposed to be geocoded, the user can define +the Area of Interest aoi in form of a GeoJSON with a single geometry. Alternatively, the +a-flag can be used to limit the geocoded area to the current computational region. This also +changes the resolution of the geocoded products to that of the current computational region.

In the nprocs option, users can specify the number of parallel process to run. If more than one diff --git a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py index 84b335e4..a3b6e82e 100755 --- a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py +++ b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py @@ -149,6 +149,16 @@ # % description: Link resulting data and do not read statistics # %end +# %flag +# % key: a +# % description: Fetch aoi and pixel spacing from current region +# %end + +# %rules +# % exclusive: aoi,-a +# %end + +import atexit import os import shutil import sys @@ -162,6 +172,19 @@ import grass.script as gs +rm_vec = [] +rm_files = [] +def cleanup(): + """Cleanup fuction (can be extended)""" + nulldev = open(os.devnull, "w", encoding="utf-8") + kwargs = {"flags": "f", "quiet": True, "stderr": nulldev} + for rmvec in rm_vec: + if gs.find_file(name=rmvec, element="vector")["file"]: + gs.run_command("g.remove", type="vector", name=rmvec, **kwargs) + for rmfile in rm_files: + if os.path.isfile(rmfile): + os.remove(rmfile) + def get_raster_gdalpath(map_name, check_linked=True): """Get get the path to a raster map that can be opened by GDAL Checks for GDAL source of linked raster data and returns those @@ -216,7 +239,7 @@ def get_aoi_geometry(geojson_file): ) ogr_layer = ogr_dataset.GetLayerByIndex(0) if ogr_layer.GetGeomType() != 3: - gs.fatal(_("GeoJSON does not contain polygons")) + gs.warning(_("GeoJSON does not contain polygons")) if ogr_layer.GetFeatureCount() > 1: gs.warning( _("GeoJSON contains more than one geometry. Using only the first one.") @@ -248,7 +271,6 @@ def get_target_geometry(s1_file, aoi=None): bpol_geom.FlattenTo2D() aoi = ogr.CreateGeometryFromWkt(aoi, s_srs) - if bpol_geom.Within(aoi): return True if not bpol_geom.Intersect(aoi): @@ -436,7 +458,7 @@ def check_files_list(file_path_list): def main(): """Do the main work""" - + global rm_vec, rm_files # Check if gpt is available if not shutil.which("gpt"): gs.fatal(_("gpt commandline tool from ESA SNAP not found on current PATH")) @@ -514,6 +536,35 @@ def main(): # Check if gpt is available gpt_options = get_gpt_options(nprocs_inner) + # use spacing from DEM or region + spacing = float(dem_info["nsres"]) + aoi = options["aoi"] if options["aoi"] else None + if flags["a"]: + region = gs.region() + if region["nsres"] != region["ewres"]: + gs.warning( + _( + "Region nsres {nsres} differs from ewres {ewres}. " + "Using {nsres}..." + ).format(nsres=region["nsres"], ewres=region["ewres"]) + ) + spacing = float(region["nsres"]) + + # save current region as geojson + region_vect = f"cur_region_{os.getpid()}" + rm_vec.append(region_vect) + gs.run_command("v.in.region", output=region_vect, quiet=True) + cur_reg_json = f"{gs.tempfile(create=False)}.geojson" + rm_files.append(cur_reg_json) + gs.run_command("v.out.ogr", input=region_vect, output=cur_reg_json, + format="GeoJSON", quiet=True) + # reproject to 4326 + reg_json_4326 = f"{gs.tempfile(create=False)}_4326.geojson" + rm_files.append(reg_json_4326) + ds = gdal.VectorTranslate(reg_json_4326, cur_reg_json, reproject=True, dstSRS="EPSG:4326") + ds = None + aoi = reg_json_4326 + # Setup function for geocoding geocode_kwargs = { "t_srs": location_crs_wkt, @@ -525,7 +576,7 @@ def main(): "speckleFilter": speckle_filter_dict[options["speckle_filter"]] if options["speckle_filter"] else None, - "spacing": float(dem_info["nsres"]), + "spacing": spacing, "externalDEMFile": str(dem_info["GDAL_path"]), "externalDEMNoDataValue": -2147483678.0 if dem_info["datatype"] == "DCELL" @@ -543,7 +594,6 @@ def main(): "gpt_args": gpt_options[1], "tmpdir": options["temporary_directory"] or tempfile.gettempdir(), } - if export_extra and "gammaSigmaRatio" in export_extra: geocode_kwargs["refarea"] = ["sigma0", "gamma0"] @@ -551,7 +601,7 @@ def main(): _geocode_snap = partial( process_image_file, kwargs=geocode_kwargs, - aoi=get_aoi_geometry(options["aoi"]) if options["aoi"] else None, + aoi=get_aoi_geometry(aoi) if aoi else None, import_flags=flags, ) @@ -602,4 +652,5 @@ def main(): ) ) + atexit.register(cleanup) sys.exit(main()) From 24a1c63215777df362e91abc814573d19ea42c7c Mon Sep 17 00:00:00 2001 From: Guido Riembauer Date: Tue, 16 Apr 2024 14:24:55 +0200 Subject: [PATCH 2/2] reactivate gs.fatal --- .../i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py index a3b6e82e..a3d1f52a 100755 --- a/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py +++ b/src/imagery/i.sentinel1.pyrosargeocode/i.sentinel1.pyrosargeocode.py @@ -239,7 +239,7 @@ def get_aoi_geometry(geojson_file): ) ogr_layer = ogr_dataset.GetLayerByIndex(0) if ogr_layer.GetGeomType() != 3: - gs.warning(_("GeoJSON does not contain polygons")) + gs.fatal(_("GeoJSON does not contain polygons")) if ogr_layer.GetFeatureCount() > 1: gs.warning( _("GeoJSON contains more than one geometry. Using only the first one.")