From 475540b5d65fadec00f1984d9669108a0485efe1 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:05:02 +0200 Subject: [PATCH 1/7] [search.combine_polygons] new function --- docs/source/api.rst | 1 + s1ard/search.py | 70 ++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 70 insertions(+), 1 deletion(-) diff --git a/docs/source/api.rst b/docs/source/api.rst index a21c987..125a19a 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -177,6 +177,7 @@ Scene Search check_acquisition_completeness collect_neighbors scene_select + combine_polygons Metadata -------- diff --git a/s1ard/search.py b/s1ard/search.py index 92c8235..bd41a8a 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -6,11 +6,12 @@ from datetime import datetime, timedelta, timezone from pystac_client import Client from pystac_client.stac_api_io import StacApiIO -from spatialist import Vector, crsConvert +from spatialist.vector import Vector, crsConvert import asf_search as asf from pyroSAR import identify_many, ID from s1ard.ancillary import buffer_time from s1ard.tile_extraction import aoi_from_tile, tile_from_aoi +from osgeo import ogr, osr class ASF(ID): @@ -698,3 +699,70 @@ def check_acquisition_completeness(archive, scenes): if len(messages) != 0: text = '\n - '.join(messages) raise RuntimeError(f'missing the following scenes:\n - {text}') + + +def combine_polygons(vector, crs=4326, multipolygon=False, layer_name='combined'): + """ + Combine polygon vector objects into one. + The output is a single vector object with the polygons either stored in + separate features or combined into a single multipolygon geometry. + + Parameters + ---------- + vector: list[spatialist.vector.Vector] + crs: int or str + the target CRS. Default: EPSG:4326 + multipolygon: bool + combine all polygons into one multipolygon? + Default False: write each polygon into a separate feature. + layer_name: str + the layer name of the output vector object. + + Returns + ------- + spatialist.vector.Vector + """ + if not isinstance(vector, list): + raise TypeError("'vectorobject' must be a list") + geometry_names = [] + for item in vector: + for feature in item.layer: + geom = feature.GetGeometryRef() + geometry_names.append(geom.GetGeometryName()) + item.layer.ResetReading() + geom = None + geometry_names = list(set(geometry_names)) + if not all(x == 'POLYGON' for x in geometry_names): + raise RuntimeError('All geometries must be of type POLYGON') + + vec = Vector(driver='Memory') + srs_out = crsConvert(crs, 'osr') + if multipolygon: + geom_type = ogr.wkbMultiPolygon + geom_out = ogr.Geometry(geom_type) + else: + geom_type = ogr.wkbPolygon + geom_out = [] + vec.addlayer(name=layer_name, srs=srs_out, geomType=geom_type) + for item in vector: + if item.srs.IsSame(srs_out): + coord_trans = None + else: + coord_trans = osr.CoordinateTransformation(item.srs, srs_out) + for feature in item.layer: + geom = feature.GetGeometryRef() + if coord_trans is not None: + geom.Transform(coord_trans) + if multipolygon: + geom_out.AddGeometry(geom.Clone()) + else: + geom_out.append(geom.Clone()) + item.layer.ResetReading() + geom = None + if multipolygon: + vec.addfeature(geom_out) + else: + for geom in geom_out: + vec.addfeature(geom) + geom_out = None + return vec From 91ccf8a603b4b622be59f272a363b482fb19d118 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:06:16 +0200 Subject: [PATCH 2/7] [search.scene_select] make use of new func combine_polygons --- s1ard/search.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/s1ard/search.py b/s1ard/search.py index bd41a8a..d2c0ebf 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -517,7 +517,6 @@ def scene_select(archive, kml_file, aoi_tiles=None, aoi_geometry=None, **kwargs) args['return_value'] = 'ASF' vec = None - selection = [] if aoi_tiles is not None: vec = aoi_from_tile(kml=kml_file, tile=aoi_tiles) elif aoi_geometry is not None: @@ -569,9 +568,8 @@ def scene_select(archive, kml_file, aoi_tiles=None, aoi_geometry=None, **kwargs) if isinstance(archive, ASFArchive): args['return_value'] = 'url' - for item in vec: - args['vectorobject'] = item - selection.extend(archive.select(**args)) + with combine_polygons(vec, multipolygon=True) as combined: + selection = archive.select(**args, vectorobject=combined) del vec, args return sorted(list(set(selection))), aoi_tiles From 858fd88441c72093d9bc8576e5d7679684928788 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:25:52 +0200 Subject: [PATCH 3/7] [search.scene_select] bug fix --- s1ard/search.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/s1ard/search.py b/s1ard/search.py index d2c0ebf..9463318 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -569,7 +569,8 @@ def scene_select(archive, kml_file, aoi_tiles=None, aoi_geometry=None, **kwargs) args['return_value'] = 'url' with combine_polygons(vec, multipolygon=True) as combined: - selection = archive.select(**args, vectorobject=combined) + args['vectorobject'] = combined + selection = archive.select(**args) del vec, args return sorted(list(set(selection))), aoi_tiles From 85e9a68147e6df568ce1abad3aac0e7c1705094c Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:46:46 +0200 Subject: [PATCH 4/7] [STACArchive.select] ensure vectorobject only has one feature --- s1ard/search.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/s1ard/search.py b/s1ard/search.py index 9463318..a1c2ba7 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -181,7 +181,7 @@ def select(self, sensor=None, product=None, acquisition_mode=None, the data take ID in decimal representation. Requires custom STAC key `s1:datatake`. vectorobject: spatialist.vector.Vector or None - a geometry with which the scenes need to overlap + a geometry with which the scenes need to overlap. The object may only contain one feature. date_strict: bool treat dates as strict limits or also allow flexible limits to incorporate scenes whose acquisition period overlaps with the defined limit? @@ -242,6 +242,8 @@ def select(self, sensor=None, product=None, acquisition_mode=None, flt['args'].append(arg) elif key == 'vectorobject': if isinstance(val, Vector): + if val.nfeatures > 1: + raise RuntimeError("'vectorobject' contains more than one feature.") with val.clone() as vec: vec.reproject(4326) ext = vec.extent From 68ba1dcb7d7cb8ae470114a19620e6c580f26ba1 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:47:13 +0200 Subject: [PATCH 5/7] [asf_select] ensure vectorobject only has one feature --- s1ard/search.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/s1ard/search.py b/s1ard/search.py index a1c2ba7..ac59a33 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -325,7 +325,7 @@ def select(sensor=None, product=None, acquisition_mode=None, mindate=None, maxdate: str or datetime.datetime or None the maximum acquisition date vectorobject: spatialist.vector.Vector or None - a geometry with which the scenes need to overlap + a geometry with which the scenes need to overlap. The object may only contain one feature. date_strict: bool treat dates as strict limits or also allow flexible limits to incorporate scenes whose acquisition period overlaps with the defined limit? @@ -368,7 +368,7 @@ def asf_select(sensor, product, acquisition_mode, mindate, maxdate, maxdate: str or datetime.datetime the maximum acquisition date vectorobject: spatialist.vector.Vector or None - a geometry with which the scenes need to overlap + a geometry with which the scenes need to overlap. The object may only contain one feature. return_value: str or list[str] the metadata return value; if `ASF`, an :class:`~s1ard.search.ASF` object is returned; further string options specify certain properties to return: `beamModeType`, `browse`, @@ -403,6 +403,8 @@ def asf_select(sensor, product, acquisition_mode, mindate, maxdate, else: beam_mode = acquisition_mode if vectorobject is not None: + if vectorobject.nfeatures > 1: + raise RuntimeError("'vectorobject' contains more than one feature.") with vectorobject.clone() as geom: geom.reproject(4326) geometry = geom.convert2wkt(set3D=False)[0] From b179359dd1fe382ea424ae8e607e8419e77f77e5 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 13:48:06 +0200 Subject: [PATCH 6/7] [STACArchive.select] search using the input geometry, not its bounding box --- s1ard/search.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/s1ard/search.py b/s1ard/search.py index ac59a33..ff476ae 100644 --- a/s1ard/search.py +++ b/s1ard/search.py @@ -246,8 +246,10 @@ def select(self, sensor=None, product=None, acquisition_mode=None, raise RuntimeError("'vectorobject' contains more than one feature.") with val.clone() as vec: vec.reproject(4326) - ext = vec.extent - args['bbox'] = [ext['xmin'], ext['ymin'], ext['xmax'], ext['ymax']] + feat = vec.getFeatureByIndex(0) + json = feat.ExportToJson(as_object=True) + feat = None + args['intersects'] = json else: raise TypeError('argument vectorobject must be of type spatialist.vector.Vector') else: From 9da66fd0fff9c46ca262aa28654b70c785efce15 Mon Sep 17 00:00:00 2001 From: John Truckenbrodt Date: Tue, 23 Jul 2024 15:39:51 +0200 Subject: [PATCH 7/7] [tile_extraction] convert all geometries to Polygon using new function multipolygon2polygon --- docs/source/api.rst | 1 + s1ard/tile_extraction.py | 45 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 42 insertions(+), 4 deletions(-) diff --git a/docs/source/api.rst b/docs/source/api.rst index 4a5b1af..05bd292 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -137,6 +137,7 @@ Tile Extraction aoi_from_scene aoi_from_tile description2dict + multipolygon2polygon tile_from_aoi Ancillary Functions diff --git a/s1ard/tile_extraction.py b/s1ard/tile_extraction.py index bdd0762..325a3fa 100644 --- a/s1ard/tile_extraction.py +++ b/s1ard/tile_extraction.py @@ -4,6 +4,7 @@ from spatialist.vector import Vector, wkt2vector, bbox from spatialist.auxil import utm_autodetect from s1ard.ancillary import get_max_ext, buffer_min_overlap, get_kml +from osgeo import ogr def tile_from_aoi(vector, epsg=None, strict=True, return_geometries=False, tilenames=None): @@ -70,15 +71,16 @@ def tile_from_aoi(vector, epsg=None, strict=True, return_geometries=False, tilen else: continue if return_geometries: + wkt = multipolygon2polygon(attrib['UTM_WKT']) if reproject: - with wkt2vector(attrib['UTM_WKT'], attrib['EPSG']) as tmp: + with wkt2vector(wkt, attrib['EPSG']) as tmp: tmp.reproject(epsg_target) ext = tmp.extent for k, v in ext.items(): ext[k] = round(v / 10) * 10 geom = bbox(ext, crs=epsg_target) else: - geom = wkt2vector(attrib['UTM_WKT'], attrib['EPSG']) + geom = wkt2vector(wkt, attrib['EPSG']) geom.mgrs = tilename tiles.append(geom) else: @@ -117,10 +119,11 @@ def aoi_from_tile(tile): feat = vec.getFeatureByAttribute('Name', tilename) attrib = description2dict(feat.GetField('Description')) feat = None + wkt = multipolygon2polygon(attrib['UTM_WKT']) if epsg is None: - return wkt2vector(attrib['UTM_WKT'], attrib['EPSG']) + return wkt2vector(wkt, attrib['EPSG']) else: - with wkt2vector(attrib['UTM_WKT'], attrib['EPSG']) as tmp: + with wkt2vector(wkt, attrib['EPSG']) as tmp: tmp.reproject(int(epsg)) ext = tmp.extent for k, v in ext.items(): @@ -230,3 +233,37 @@ def fn(x): out.append({'extent': ext, 'epsg': epsg, 'align_x': align_x, 'align_y': align_y}) return out + + +def multipolygon2polygon(wkt): + """ + Convert a MultiPolygon WKT with one Polygon to a simple Polygon WKT. + The Sentinel-2 KML grid file stores all geometries as MultiPolygons. + This function simply converts the geometries to simple Polygons. + Not all geometries in the KML file have been checked. In case there are + ever multiple Polygons in one MultiPolygon, an RuntimeError is raised. + All other geometries are returned as is. + + Parameters + ---------- + wkt: str + A geometry WKT representation. + + Returns + ------- + str + the output WKT geometry. Either the original geometry or + a Polygon extracted from a single-Polygon MultiPolygon. + """ + geom1 = ogr.CreateGeometryFromWkt(wkt) + poly_type = geom1.GetGeometryType() + poly_num = geom1.GetGeometryCount() + if poly_type == ogr.wkbMultiPolygon: + if poly_num == 1: + geom2 = geom1.GetGeometryRef(0).Clone() + wkt = geom2.ExportToWkt() + geom2 = None + else: + raise RuntimeError('got a MultiPolygon with multiple Polygons') + geom1 = None + return wkt