diff --git a/docs/source/api.rst b/docs/source/api.rst index c7339f1..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 @@ -180,6 +181,7 @@ Scene Search check_acquisition_completeness collect_neighbors scene_select + combine_polygons Metadata -------- diff --git a/s1ard/search.py b/s1ard/search.py index c69863f..6d73a04 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): @@ -180,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? @@ -241,10 +242,14 @@ 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 - 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: @@ -322,7 +327,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? @@ -365,7 +370,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`, @@ -400,6 +405,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] @@ -514,7 +521,6 @@ def scene_select(archive, aoi_tiles=None, aoi_geometry=None, **kwargs): args['return_value'] = 'ASF' vec = None - selection = [] if aoi_tiles is not None: vec = aoi_from_tile(tile=aoi_tiles) elif aoi_geometry is not None: @@ -566,9 +572,9 @@ def scene_select(archive, 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: + args['vectorobject'] = combined + selection = archive.select(**args) del vec, args return sorted(list(set(selection))), aoi_tiles @@ -696,3 +702,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 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