diff --git a/dep_tools/grids.py b/dep_tools/grids.py index 73e8ec2..ec5c3a8 100644 --- a/dep_tools/grids.py +++ b/dep_tools/grids.py @@ -115,7 +115,6 @@ def grid( Otherwise, a GeoDataFrame containing only the portions of each tile that intersect the given GeoDataFrame is returned. """ - if intersect_with is not None: if return_type != "GridSpec": full_grid = _geoseries(resolution, crs) diff --git a/dep_tools/landsat_utils.py b/dep_tools/landsat_utils.py index 74d35f1..f0b8ae3 100644 --- a/dep_tools/landsat_utils.py +++ b/dep_tools/landsat_utils.py @@ -6,7 +6,7 @@ from shapely.geometry import box from xarray import DataArray, Dataset -from dep_tools.utils import bbox_across_180 +from dep_tools.utils import bbox_across_180, fix_winding def cloud_mask( @@ -66,13 +66,20 @@ def mask_clouds( return xr.where(~mask) +def _pathrows(): + pathrows = GeoDataFrame( + read_file( + "https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/WRS2_descending_0.zip" + ) + ) + pathrows["geometry"] = pathrows.geometry.apply(fix_winding) + return pathrows + + def pathrows_in_area(area: GeoDataFrame, pathrows: GeoDataFrame | None = None): if pathrows is None: - pathrows = GeoDataFrame( - read_file( - "https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/WRS2_descending_0.zip" - ) - ) + pathrows = _pathrows() + bbox = bbox_across_180(area) if isinstance(bbox, tuple): return pathrows[ @@ -99,9 +106,7 @@ def items_in_pathrows( def pathrow_with_greatest_area(shapes: GeoDataFrame) -> Tuple[str, str]: - pathrows = read_file( - "https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/WRS2_descending_0.zip" - ) + pathrows = _pathrows() intersection = shapes.overlay(pathrows, how="intersection") row_with_greatest_area = intersection.iloc[[intersection.geometry.area.idxmax()]] return (row_with_greatest_area.PATH.item(), row_with_greatest_area.ROW.item()) diff --git a/dep_tools/searchers.py b/dep_tools/searchers.py index e6f00b5..38406a6 100644 --- a/dep_tools/searchers.py +++ b/dep_tools/searchers.py @@ -34,7 +34,9 @@ class PystacSearcher(Searcher): :func:`dep_tools.utils.search_across_180`. Args: - client: A search client. + catalog: The URL of a stac catalog. if client is specified, this is + ignored. + client: A search client. Either this or catalog must be specified. raise_empty_collection_error: Whether an EmptyCollectionError exception should be returned if no stac items are found. **kwargs: Additional arguments passed to client.search(). For example, @@ -44,11 +46,20 @@ class PystacSearcher(Searcher): def __init__( self, - catalog: str, + catalog: str | None = None, + client: Client | None = None, raise_empty_collection_error: bool = True, **kwargs, ): - self._client = Client.open(catalog) + if client and catalog: + warnings.warn( + "Arguments for both 'client' and 'catalog' passed to PystacSearcher, ignoring catalog" + ) + + if not (client or catalog): + raise ValueError("Must specify either client or catalog") + + self._client = client if client else Client.open(catalog) self._raise_errors = raise_empty_collection_error self._kwargs = kwargs @@ -65,12 +76,12 @@ def search(self, area: GeoDataFrame | GeoBox) -> ItemCollection: region=area, client=self._client, **self._kwargs ) - if len(item_collection) == 0 and self._raise_errors: - raise EmptyCollectionError() - fix_bad_epsgs(item_collection) item_collection = remove_bad_items(item_collection) + if len(item_collection) == 0 and self._raise_errors: + raise EmptyCollectionError() + return item_collection @@ -80,7 +91,9 @@ class LandsatPystacSearcher(PystacSearcher): query, just use :class:PystacSearcher. Args: - client: A search client. + catalog: The URL of a stac catalog. if client is specified, this is + ignored. + client: A search client. Either this or catalog must be specified. raise_empty_collection_error: Whether an EmptyCollectionError exception should be returned if no stac items are found. search_intersecting_pathrows: Whether to use landsat pathrows which @@ -99,7 +112,8 @@ class LandsatPystacSearcher(PystacSearcher): def __init__( self, - catalog: str = "https://planetarycomputer.microsoft.com/api/stac/v1/", + catalog: str | None = None, + client: Client | None = None, collections: list[str] | None = ["landsat-c2-l2"], raise_empty_collection_error: bool = True, search_intersecting_pathrows: bool = False, @@ -110,6 +124,7 @@ def __init__( ): super().__init__( catalog=catalog, + client=client, raise_empty_collection_error=raise_empty_collection_error, **kwargs, ) @@ -140,16 +155,9 @@ def __init__( self._kwargs["query"] = query - if self._search_intersecting_pathrows: - self._landsat_pathrows = read_file( - "https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/WRS2_descending_0.zip" - ) - def search(self, area: GeoDataFrame): search_area = ( - pathrows_in_area(area, self._landsat_pathrows) - if self._search_intersecting_pathrows - else area + pathrows_in_area(area) if self._search_intersecting_pathrows else area ) try: items = super().search(search_area) @@ -164,4 +172,7 @@ def search(self, area: GeoDataFrame): if self._search_intersecting_pathrows: items = items_in_pathrows(items, search_area) + if len(items) == 0 and self._raise_errors: + raise EmptyCollectionError() + return items diff --git a/dep_tools/utils.py b/dep_tools/utils.py index cb7e70b..43be0fe 100644 --- a/dep_tools/utils.py +++ b/dep_tools/utils.py @@ -17,7 +17,13 @@ from pystac import ItemCollection import pystac_client from retry import retry -from shapely.geometry import LineString, MultiLineString, GeometryCollection +from shapely.geometry import ( + LineString, + MultiLineString, + GeometryCollection, + MultiPolygon, +) +from shapely.geometry.polygon import orient from xarray import DataArray, Dataset # Set the timeout to five minutes, which is an extremely long time @@ -108,6 +114,23 @@ def bbox_across_180(region: GeoDataFrame | GeoBox) -> BBOX | tuple[BBOX, BBOX]: return BBOX(bbox) +from shapely import Geometry + + +def fix_winding(geom: Geometry) -> Geometry: + """Fixes the orientation of the exterior coordinates of Polygon and + MultiPolygon geometry, which should run counterclockwise. + """ + # This is non-essential but resolves a barrage of warnings from + # the antimeridian package + if geom.geom_type == "Polygon" and not geom.exterior.is_ccw: + return orient(geom) + elif geom.geom_type == "MultiPolygon": + return MultiPolygon([fix_winding(p) for p in geom.geoms]) + else: + return geom + + def _fix_geometry(geometry): if isinstance(geometry, GeometryCollection): return GeometryCollection( diff --git a/tests/test_fix_winding.py b/tests/test_fix_winding.py new file mode 100644 index 0000000..7630db4 --- /dev/null +++ b/tests/test_fix_winding.py @@ -0,0 +1,8 @@ +from shapely.geometry import Polygon, MultiPolygon +from dep_tools.utils import fix_winding + + +def test_fix_winding(): + bad_poly = Polygon(shell=[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]) + assert not bad_poly.exterior.is_ccw + assert fix_winding(bad_poly).exterior.is_ccw diff --git a/tests/test_landsat_mspc_stac_api.py b/tests/test_landsat_mspc_stac_api.py index f4e4e38..7579433 100644 --- a/tests/test_landsat_mspc_stac_api.py +++ b/tests/test_landsat_mspc_stac_api.py @@ -29,6 +29,7 @@ def test_hobart_all(hobart_area): exclude_platforms=None, only_tier_one=False, fall_back_to_tier_two=False, + catalog="https://planetarycomputer.microsoft.com/api/stac/v1/", ) items = list(searcher.search(hobart_area)) @@ -42,6 +43,7 @@ def test_hobart_tier_one_only(hobart_area): exclude_platforms=None, only_tier_one=True, fall_back_to_tier_two=False, + catalog="https://planetarycomputer.microsoft.com/api/stac/v1/", ) items = list(searcher.search(hobart_area)) @@ -55,6 +57,7 @@ def test_hobart_tier_exclude_7(hobart_area): exclude_platforms=["landsat-7"], only_tier_one=False, fall_back_to_tier_two=False, + catalog="https://planetarycomputer.microsoft.com/api/stac/v1/", ) items = list(searcher.search(hobart_area)) diff --git a/tests/test_landsat_mspc_stac_search_fixes.py b/tests/test_landsat_mspc_stac_search_fixes.py index e891dad..111423e 100644 --- a/tests/test_landsat_mspc_stac_search_fixes.py +++ b/tests/test_landsat_mspc_stac_search_fixes.py @@ -27,7 +27,9 @@ def test_pathrows_in_area(near_antimeridian_area): def test_search_for_stac_items_with_bad_geoms(near_antimeridian_area): searcher = LandsatPystacSearcher( - datetime=DATETIME, search_intersecting_pathrows=True + datetime=DATETIME, + search_intersecting_pathrows=True, + catalog="https://planetarycomputer.microsoft.com/api/stac/v1/", ) items = searcher.search(near_antimeridian_area) diff --git a/tests/test_searchers.py b/tests/test_searchers.py index b0f347c..f27fddd 100644 --- a/tests/test_searchers.py +++ b/tests/test_searchers.py @@ -36,7 +36,10 @@ def test_PystacSearcher(area, mspc_catalog): def test_LandsatPystacSearcher_exclude_platforms(area): - s = LandsatPystacSearcher(exclude_platforms=["landsat-7"]) + s = LandsatPystacSearcher( + exclude_platforms=["landsat-7"], + catalog="https://planetarycomputer.microsoft.com/api/stac/v1/", + ) s.search(area)