Skip to content

Commit

Permalink
Merge branch 'winding-fix' into retry-search
Browse files Browse the repository at this point in the history
  • Loading branch information
jessjaco committed Oct 25, 2024
2 parents d7eae37 + 587560c commit b850551
Show file tree
Hide file tree
Showing 5 changed files with 53 additions and 22 deletions.
1 change: 0 additions & 1 deletion dep_tools/grids.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
23 changes: 14 additions & 9 deletions dep_tools/landsat_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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[
Expand All @@ -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())
18 changes: 7 additions & 11 deletions dep_tools/searchers.py
Original file line number Diff line number Diff line change
Expand Up @@ -76,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


Expand Down Expand Up @@ -155,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)
Expand All @@ -179,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
25 changes: 24 additions & 1 deletion dep_tools/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(
Expand Down
8 changes: 8 additions & 0 deletions tests/test_fix_winding.py
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit b850551

Please sign in to comment.