Skip to content

Commit

Permalink
Merge pull request #214 from SAR-ARD/feature/search_multipolygon
Browse files Browse the repository at this point in the history 
improve geometry handling in search queries
  • Loading branch information
@johntruckenbrodt
johntruckenbrodt authored Sep 9, 2024
2 parents 992e267 + 9da66fd commit 0abc5ec
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 14 deletions.
2 changes: 2 additions & 0 deletions docs/source/api.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,7 @@ Tile Extraction
aoi_from_scene
aoi_from_tile
description2dict
multipolygon2polygon
tile_from_aoi

Ancillary Functions
Expand Down Expand Up @@ -180,6 +181,7 @@ Scene Search
check_acquisition_completeness
collect_neighbors
scene_select
combine_polygons

Metadata
--------
Expand Down
93 changes: 83 additions & 10 deletions s1ard/search.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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`,
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
45 changes: 41 additions & 4 deletions s1ard/tile_extraction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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():
Expand Down Expand Up @@ -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

0 comments on commit 0abc5ec

Please sign in to comment.