|
| 1 | +import abc |
| 2 | +import math |
| 3 | +from typing import Dict, List, NamedTuple, Optional, Union |
| 4 | + |
| 5 | +import shapely |
| 6 | +from shapely.geometry import MultiPolygon, Polygon |
| 7 | + |
| 8 | +from openeo.util import normalize_crs |
| 9 | + |
| 10 | + |
| 11 | +class JobSplittingFailure(Exception): |
| 12 | + pass |
| 13 | + |
| 14 | + |
| 15 | +class _BoundingBox(NamedTuple): |
| 16 | + """Simple NamedTuple container for a bounding box""" |
| 17 | + |
| 18 | + # TODO: this should be moved to more general utility module, and/or merged with existing BBoxDict |
| 19 | + |
| 20 | + west: float |
| 21 | + south: float |
| 22 | + east: float |
| 23 | + north: float |
| 24 | + crs: int = 4326 |
| 25 | + |
| 26 | + @classmethod |
| 27 | + def from_dict(cls, d: Dict) -> "_BoundingBox": |
| 28 | + """Create a bounding box from a dictionary""" |
| 29 | + if d.get("crs") is not None: |
| 30 | + d["crs"] = normalize_crs(d["crs"]) |
| 31 | + return cls(**{k: d[k] for k in cls._fields if k not in cls._field_defaults or k in d}) |
| 32 | + |
| 33 | + @classmethod |
| 34 | + def from_polygon(cls, polygon: Union[MultiPolygon, Polygon], crs: Optional[int] = None) -> "_BoundingBox": |
| 35 | + """Create a bounding box from a shapely Polygon or MultiPolygon""" |
| 36 | + crs = normalize_crs(crs) |
| 37 | + return cls(*polygon.bounds, crs=4326 if crs is None else crs) |
| 38 | + |
| 39 | + def as_dict(self) -> Dict: |
| 40 | + return self._asdict() |
| 41 | + |
| 42 | + def as_polygon(self) -> Polygon: |
| 43 | + """Get bounding box as a shapely Polygon""" |
| 44 | + return shapely.geometry.box(minx=self.west, miny=self.south, maxx=self.east, maxy=self.north) |
| 45 | + |
| 46 | + |
| 47 | +class _TileGridInterface(metaclass=abc.ABCMeta): |
| 48 | + """Interface for tile grid classes""" |
| 49 | + |
| 50 | + @abc.abstractmethod |
| 51 | + # TODO: is it intentional that this method returns a list of non-multi polygons even if the input can be multi-polygon? |
| 52 | + # TODO: typehint states that geometry can be a dict too, but that is very liberal, it's probably just about bounding box kind of dicts? |
| 53 | + def get_tiles(self, geometry: Union[Dict, MultiPolygon, Polygon]) -> List[Polygon]: |
| 54 | + """Calculate tiles to cover given bounding box""" |
| 55 | + ... |
| 56 | + |
| 57 | + |
| 58 | +class _SizeBasedTileGrid(_TileGridInterface): |
| 59 | + """ |
| 60 | + Specification of a tile grid, parsed from a size and a projection. |
| 61 | + The size is in m for UTM projections or degrees for WGS84. |
| 62 | + """ |
| 63 | + |
| 64 | + def __init__(self, *, epsg: int, size: float): |
| 65 | + # TODO: normalize_crs does not necessarily return an int (could also be a WKT2 string, or even None), but further logic seems to assume it's an int |
| 66 | + self.epsg = normalize_crs(epsg) |
| 67 | + self.size = size |
| 68 | + |
| 69 | + @classmethod |
| 70 | + def from_size_projection(cls, *, size: float, projection: str) -> "_SizeBasedTileGrid": |
| 71 | + """Create a tile grid from size and projection""" |
| 72 | + # TODO: the constructor also does normalize_crs, so this factory looks like overkill at the moment |
| 73 | + return cls(epsg=normalize_crs(projection), size=size) |
| 74 | + |
| 75 | + def _epsg_is_meters(self) -> bool: |
| 76 | + """Check if the projection unit is in meters. (EPSG:3857 or UTM)""" |
| 77 | + # TODO: this is a bit misleading: this code just checks some EPSG ranges (UTM and 3857) and calls all the rest to be not in meters. |
| 78 | + # It would be better to raise an exception on unknown EPSG codes than claiming they're not in meter |
| 79 | + return 32601 <= self.epsg <= 32660 or 32701 <= self.epsg <= 32760 or self.epsg == 3857 |
| 80 | + |
| 81 | + @staticmethod |
| 82 | + def _split_bounding_box(to_cover: _BoundingBox, x_offset: float, tile_size: float) -> List[Polygon]: |
| 83 | + """ |
| 84 | + Split a bounding box into tiles of given size and projection. |
| 85 | + :param to_cover: bounding box dict with keys "west", "south", "east", "north", "crs" |
| 86 | + :param x_offset: offset to apply to the west and east coordinates |
| 87 | + :param tile_size: size of tiles in unit of measure of the projection |
| 88 | + :return: list of tiles (polygons) |
| 89 | + """ |
| 90 | + xmin = int(math.floor((to_cover.west - x_offset) / tile_size)) |
| 91 | + xmax = int(math.ceil((to_cover.east - x_offset) / tile_size)) - 1 |
| 92 | + ymin = int(math.floor(to_cover.south / tile_size)) |
| 93 | + ymax = int(math.ceil(to_cover.north / tile_size)) - 1 |
| 94 | + |
| 95 | + tiles = [] |
| 96 | + for x in range(xmin, xmax + 1): |
| 97 | + for y in range(ymin, ymax + 1): |
| 98 | + tiles.append( |
| 99 | + _BoundingBox( |
| 100 | + west=max(x * tile_size + x_offset, to_cover.west), |
| 101 | + south=max(y * tile_size, to_cover.south), |
| 102 | + east=min((x + 1) * tile_size + x_offset, to_cover.east), |
| 103 | + north=min((y + 1) * tile_size, to_cover.north), |
| 104 | + ).as_polygon() |
| 105 | + ) |
| 106 | + |
| 107 | + return tiles |
| 108 | + |
| 109 | + def get_tiles(self, geometry: Union[Dict, MultiPolygon, Polygon]) -> List[Polygon]: |
| 110 | + if isinstance(geometry, dict): |
| 111 | + bbox = _BoundingBox.from_dict(geometry) |
| 112 | + |
| 113 | + elif isinstance(geometry, Polygon) or isinstance(geometry, MultiPolygon): |
| 114 | + bbox = _BoundingBox.from_polygon(geometry, crs=self.epsg) |
| 115 | + |
| 116 | + else: |
| 117 | + raise JobSplittingFailure("geometry must be a dict or a shapely.geometry.Polygon or MultiPolygon") |
| 118 | + |
| 119 | + # TODO: being a meter based EPSG does not imply that offset should be 500_000 |
| 120 | + x_offset = 500_000 if self._epsg_is_meters() else 0 |
| 121 | + |
| 122 | + tiles = _SizeBasedTileGrid._split_bounding_box(to_cover=bbox, x_offset=x_offset, tile_size=self.size) |
| 123 | + |
| 124 | + return tiles |
| 125 | + |
| 126 | + |
| 127 | +def split_area( |
| 128 | + aoi: Union[Dict, MultiPolygon, Polygon], projection: str = "EPSG:3857", tile_size: float = 20_000.0 |
| 129 | +) -> List[Polygon]: |
| 130 | + """ |
| 131 | + Split area of interest into tiles of given size and projection. |
| 132 | + :param aoi: area of interest (bounding box or shapely polygon) |
| 133 | + :param projection: projection to use for splitting. Default is web mercator (EPSG:3857) |
| 134 | + :param tile_size: size of tiles in unit of measure of the projection |
| 135 | + :return: list of tiles (polygons). |
| 136 | + """ |
| 137 | + # TODO EPSG 3857 is probably not a good default projection. Probably better to make it a required parameter |
| 138 | + if isinstance(aoi, dict): |
| 139 | + # TODO: this possibly overwrites the given projection without the user noticing, making usage confusing |
| 140 | + projection = aoi.get("crs", projection) |
| 141 | + |
| 142 | + tile_grid = _SizeBasedTileGrid.from_size_projection(size=tile_size, projection=projection) |
| 143 | + return tile_grid.get_tiles(aoi) |
0 commit comments