Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,10 @@ <h2>DESCRIPTION</h2>
the <b>l-flag</b>, <b>m-flag</b>, or <b>r-flag</b> are given.

<p>
If only a specifc extent of the Sentinel-1 images is supposed to be geocoded, the user can define
the Area of Interest <b>aoi</b> 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 <b>aoi</b> in form of a GeoJSON with a single geometry. Alternatively, the
<b>a-flag</b> 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.

<p>
In the <b>nprocs</b> option, users can specify the number of parallel process to run. If more than one
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
import os
import json
import os

import shutil
import sys
Expand All @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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"))
Expand Down Expand Up @@ -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
Comment on lines +554 to +566
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you consider using g.region -b to get region bounds in lat/long directly and and then write it as WKT to GeoJSON? That would mean less temporary files and maps...

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! I considered this too but the suggested way in this PR seemed less "invasive" as a WKT generated from region -b directly would need to be passed to get_target_geometry, which would mean passing more optional stuff through functions. The way it is implemented in this PR, the new region AOI geojson is created once in the beginning and all the remaining pipeline remains untouched.
So I suppose there is no strong reason to do it this way, but just the attempt to alter the existing code as little as possible.

Comment on lines +554 to +566
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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
reg_bounds = gs.parse_command("g.region", flags="lg")
reg_json_4326 = Path(f"{gs.tempfile(create=False)}_4326.geojson")
reg_json_4326.write_text(
json.dumps(
{
"type": "FeatureCollection",
"name": "region_aoi",
"crs": {
"type": "name",
"properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"},
},
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
[
reg_bounds["nw_long"],
reg_bounds["nw_lat"],
],
[
reg_bounds["sw_long"],
reg_bounds["sw_lat"],
],
[
reg_bounds["se_long"],
reg_bounds["se_lat"],
],
[
reg_bounds["ne_long"],
reg_bounds["ne_lat"],
],
[
reg_bounds["nw_long"],
reg_bounds["nw_lat"],
],
]
]
],
},
}
],
},
indent=2,
),
encoding="UTF8",
)
aoi = str(reg_json_4326)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here is my line of thinking. Do you think this may work?
Then, atexit procedure might not be needed either...


# Setup function for geocoding
geocode_kwargs = {
"t_srs": location_crs_wkt,
Expand All @@ -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"
Expand All @@ -543,15 +594,14 @@ 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"]

# Pre-configure geocode function
_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,
)

Expand Down Expand Up @@ -602,4 +652,5 @@ def main():
)
)

atexit.register(cleanup)
sys.exit(main())