1""" 2Module for working with large geographical areas 3""" 4import os 5import itertools 6from abc import ABC, abstractmethod 7import json 8import math 9 10import shapely.ops 11import shapely.geometry 12from shapely.geometry import Polygon, MultiPolygon, GeometryCollection 13 14from .config import SHConfig 15from .constants import CRS 16from .data_collections import DataCollection, handle_deprecated_data_source 17from .geometry import BBox, BBoxCollection, BaseGeometry, Geometry 18from .geo_utils import transform_point 19from .ogc import WebFeatureService 20from .sentinelhub_batch import SentinelHubBatch 21 22 23class AreaSplitter(ABC): 24 """ Abstract class for splitter classes. It implements common methods used for splitting large area into smaller 25 parts. 26 """ 27 def __init__(self, shape_list, crs, reduce_bbox_sizes=False): 28 """ 29 :param shape_list: A list of geometrical shapes describing the area of interest 30 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 31 :param crs: Coordinate reference system of the shapes in `shape_list` 32 :type crs: CRS 33 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly fit 34 the given geometry in `shape_list`. 35 :type reduce_bbox_sizes: bool 36 """ 37 self.crs = CRS(crs) 38 self._parse_shape_list(shape_list, self.crs) 39 self.shape_list = shape_list 40 self.area_shape = self._join_shape_list(shape_list) 41 self.reduce_bbox_sizes = reduce_bbox_sizes 42 43 self.area_bbox = self.get_area_bbox() 44 self.bbox_list = None 45 self.info_list = None 46 47 @staticmethod 48 def _parse_shape_list(shape_list, crs): 49 """ Checks if the given list of shapes is in correct format and parses geometry objects 50 51 :param shape_list: The parameter `shape_list` from class initialization 52 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 53 :raises: ValueError 54 """ 55 if not isinstance(shape_list, list): 56 raise ValueError('Splitter must be initialized with a list of shapes') 57 58 return [AreaSplitter._parse_shape(shape, crs) for shape in shape_list] 59 60 @staticmethod 61 def _parse_shape(shape, crs): 62 """ Helper method for parsing input shapes 63 """ 64 if isinstance(shape, (Polygon, MultiPolygon)): 65 return shape 66 if isinstance(shape, BaseGeometry): 67 return shape.transform(crs).geometry 68 raise ValueError(f'The list of shapes must contain shapes of types {Polygon}, {MultiPolygon} or subtype of ' 69 f'{BaseGeometry}') 70 71 @staticmethod 72 def _parse_split_parameters(split_parameter, allow_float=False): 73 """ Parses the parameters defining the splitting of the BBox 74 75 :param split_parameter: The parameters defining the split. A tuple of int for `BBoxSplitter`, a tuple of float 76 for `BaseUtmSplitter` 77 :type split_parameter: int or (int, int) or float or (float, float) 78 :param allow_float: Whether to check for floats or not 79 :type allow_float: bool 80 :return: A tuple of n 81 :rtype: (int, int) 82 :raises: ValueError 83 """ 84 parameters_type = (int, float) if allow_float else int 85 if isinstance(split_parameter, parameters_type): 86 return split_parameter, split_parameter 87 88 if isinstance(split_parameter, (tuple, list)) and len(split_parameter) == 2 and \ 89 all(isinstance(param, parameters_type) for param in split_parameter): 90 return split_parameter[0], split_parameter[1] 91 92 extra_type = '/float' if allow_float else '' 93 raise ValueError(f'Split parameter must be an int{extra_type} or a tuple of 2 int{extra_type} but ' 94 f'{split_parameter} was given') 95 96 @staticmethod 97 def _join_shape_list(shape_list): 98 """ Joins a list of shapes together into one shape 99 100 :param shape_list: A list of geometrical shapes describing the area of interest 101 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 102 :return: A multipolygon which is a union of shapes in given list 103 :rtype: shapely.geometry.multipolygon.MultiPolygon 104 """ 105 return shapely.ops.cascaded_union(shape_list) 106 107 @abstractmethod 108 def _make_split(self): 109 """ The abstract method where the splitting will happen 110 """ 111 raise NotImplementedError 112 113 def get_bbox_list(self, crs=None, buffer=None, reduce_bbox_sizes=None): 114 """ Returns a list of bounding boxes that are the result of the split 115 116 :param crs: Coordinate reference system in which the bounding boxes should be returned. If `None` the CRS will 117 be the default CRS of the splitter. 118 :type crs: CRS or None 119 :param buffer: A percentage of each BBox size increase. This will cause neighbouring bounding boxes to overlap. 120 :type buffer: float or None 121 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly 122 fit the given geometry in `shape_list`. This overrides the same parameter from constructor 123 :type reduce_bbox_sizes: bool 124 :return: List of bounding boxes 125 :rtype: list(BBox) 126 """ 127 bbox_list = self.bbox_list 128 if buffer: 129 bbox_list = [bbox.buffer(buffer) for bbox in bbox_list] 130 131 if reduce_bbox_sizes is None: 132 reduce_bbox_sizes = self.reduce_bbox_sizes 133 if reduce_bbox_sizes: 134 bbox_list = self._reduce_sizes(bbox_list) 135 136 if crs: 137 return [bbox.transform(crs) for bbox in bbox_list] 138 return bbox_list 139 140 def get_geometry_list(self): 141 """ For each bounding box an intersection with the shape of entire given area is calculated. CRS of the returned 142 shapes is the same as CRS of the given area. 143 144 :return: List of polygons or multipolygons corresponding to the order of bounding boxes 145 :rtype: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 146 """ 147 return [self._intersection_area(bbox) for bbox in self.bbox_list] 148 149 def get_info_list(self): 150 """ Returns a list of dictionaries containing information about bounding boxes obtained in split. The order in 151 the list matches the order of the list of bounding boxes. 152 153 :return: List of dictionaries 154 :rtype: list(BBox) 155 """ 156 return self.info_list 157 158 def get_area_shape(self): 159 """ Returns a single shape of entire area described with `shape_list` parameter 160 161 :return: A multipolygon which is a union of shapes describing the area 162 :rtype: shapely.geometry.multipolygon.MultiPolygon 163 """ 164 return self.area_shape 165 166 def get_area_bbox(self, crs=None): 167 """ Returns a bounding box of the entire area 168 169 :param crs: Coordinate reference system in which the bounding box should be returned. If `None` the CRS will 170 be the default CRS of the splitter. 171 :type crs: CRS or None 172 :return: A bounding box of the area defined by the `shape_list` 173 :rtype: BBox 174 """ 175 bbox_list = [BBox(shape.bounds, crs=self.crs) for shape in self.shape_list] 176 area_minx = min([bbox.lower_left[0] for bbox in bbox_list]) 177 area_miny = min([bbox.lower_left[1] for bbox in bbox_list]) 178 area_maxx = max([bbox.upper_right[0] for bbox in bbox_list]) 179 area_maxy = max([bbox.upper_right[1] for bbox in bbox_list]) 180 bbox = BBox([area_minx, area_miny, area_maxx, area_maxy], crs=self.crs) 181 if crs is None: 182 return bbox 183 return bbox.transform(crs) 184 185 def _intersects_area(self, bbox): 186 """ Checks if the bounding box intersects the entire area 187 188 :param bbox: A bounding box 189 :type bbox: BBox 190 :return: `True` if bbox intersects the entire area else False 191 :rtype: bool 192 """ 193 return self._bbox_to_area_polygon(bbox).intersects(self.area_shape) 194 195 def _intersection_area(self, bbox): 196 """ Calculates the intersection of a given bounding box and the entire area 197 198 :param bbox: A bounding box 199 :type bbox: BBox 200 :return: A shape of intersection 201 :rtype: shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon 202 """ 203 return self._bbox_to_area_polygon(bbox).intersection(self.area_shape) 204 205 def _bbox_to_area_polygon(self, bbox): 206 """ Transforms bounding box into a polygon object in the area CRS. 207 208 :param bbox: A bounding box 209 :type bbox: BBox 210 :return: A polygon 211 :rtype: shapely.geometry.polygon.Polygon 212 """ 213 projected_bbox = bbox.transform(self.crs) 214 return projected_bbox.geometry 215 216 def _reduce_sizes(self, bbox_list): 217 """ Reduces sizes of bounding boxes 218 """ 219 return [BBox(self._intersection_area(bbox).bounds, self.crs).transform(bbox.crs) for bbox in bbox_list] 220 221 222class BBoxSplitter(AreaSplitter): 223 """ A tool that splits the given area into smaller parts. Given the area it calculates its bounding box and splits 224 it into smaller bounding boxes of equal size. Then it filters out the bounding boxes that do not intersect the 225 area. If specified by user it can also reduce the sizes of the remaining bounding boxes to best fit the area. 226 """ 227 def __init__(self, shape_list, crs, split_shape, **kwargs): 228 """ 229 :param shape_list: A list of geometrical shapes describing the area of interest 230 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 231 :param crs: Coordinate reference system of the shapes in `shape_list` 232 :type crs: CRS 233 :param split_shape: Parameter that describes the shape in which the area bounding box will be split. 234 It can be a tuple of the form `(n, m)` which means the area bounding box will be split into `n` columns 235 and `m` rows. It can also be a single integer `n` which is the same as `(n, n)`. 236 :type split_shape: int or (int, int) 237 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly fit 238 the given area geometry from `shape_list`. 239 :type reduce_bbox_sizes: bool 240 """ 241 super().__init__(shape_list, crs, **kwargs) 242 243 self.split_shape = self._parse_split_parameters(split_shape) 244 245 self._make_split() 246 247 def _make_split(self): 248 """ This method makes the split 249 """ 250 columns, rows = self.split_shape 251 bbox_partition = self.area_bbox.get_partition(num_x=columns, num_y=rows) 252 253 self.bbox_list = [] 254 self.info_list = [] 255 for i, j in itertools.product(range(columns), range(rows)): 256 if self._intersects_area(bbox_partition[i][j]): 257 self.bbox_list.append(bbox_partition[i][j]) 258 259 info = {'parent_bbox': self.area_bbox, 260 'index_x': i, 261 'index_y': j} 262 self.info_list.append(info) 263 264 265class OsmSplitter(AreaSplitter): 266 """ A tool that splits the given area into smaller parts. For the splitting it uses Open Street Map (OSM) grid on 267 the specified zoom level. It calculates bounding boxes of all OSM tiles that intersect the area. If specified by 268 user it can also reduce the sizes of the remaining bounding boxes to best fit the area. 269 """ 270 _POP_WEB_MAX = None 271 272 def __init__(self, shape_list, crs, zoom_level, **kwargs): 273 """ 274 :param shape_list: A list of geometrical shapes describing the area of interest 275 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 276 :param crs: Coordinate reference system of the shapes in `shape_list` 277 :type crs: CRS 278 :param zoom_level: A zoom level defined by OSM. Level 0 is entire world, level 1 splits the world into 279 4 parts, etc. 280 :type zoom_level: int 281 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly fit 282 the given area geometry from `shape_list`. 283 :type reduce_bbox_sizes: bool 284 """ 285 if self._POP_WEB_MAX is None: 286 OsmSplitter._POP_WEB_MAX = transform_point((180, 0), CRS.WGS84, CRS.POP_WEB)[0] 287 288 super().__init__(shape_list, crs, **kwargs) 289 self.zoom_level = zoom_level 290 291 self._make_split() 292 293 def _make_split(self, ): 294 """This method makes the split 295 """ 296 self.area_bbox = self.get_area_bbox(CRS.POP_WEB) 297 self._check_area_bbox() 298 299 self.bbox_list = [] 300 self.info_list = [] 301 self._recursive_split(self.get_world_bbox(), 0, 0, 0) 302 303 for i, bbox in enumerate(self.bbox_list): 304 self.bbox_list[i] = bbox.transform(self.crs) 305 306 def _check_area_bbox(self): 307 """ The method checks if the area bounding box is completely inside the OSM grid. That means that its latitudes 308 must be contained in the interval (-85.0511, 85.0511) 309 310 :raises: ValueError 311 """ 312 for coord in self.area_bbox: 313 if abs(coord) > self._POP_WEB_MAX: 314 raise ValueError('OsmTileSplitter only works for areas which have latitude in interval ' 315 '(-85.0511, 85.0511)') 316 317 def get_world_bbox(self): 318 """ Creates a bounding box of the entire world in EPSG: 3857 319 320 :return: Bounding box of entire world 321 :rtype: BBox 322 """ 323 return BBox((-self._POP_WEB_MAX, -self._POP_WEB_MAX, self._POP_WEB_MAX, self._POP_WEB_MAX), crs=CRS.POP_WEB) 324 325 def _recursive_split(self, bbox, zoom_level, column, row): 326 """ Method that recursively creates bounding boxes of OSM grid that intersect the area. 327 328 :param bbox: Bounding box 329 :type bbox: BBox 330 :param zoom_level: OSM zoom level 331 :type zoom_level: int 332 :param column: Column in the OSM grid 333 :type column: int 334 :param row: Row in the OSM grid 335 :type row: int 336 """ 337 if zoom_level == self.zoom_level: 338 self.bbox_list.append(bbox) 339 self.info_list.append({'zoom_level': zoom_level, 340 'index_x': column, 341 'index_y': row}) 342 return 343 344 bbox_partition = bbox.get_partition(num_x=2, num_y=2) 345 for i, j in itertools.product(range(2), range(2)): 346 if self._intersects_area(bbox_partition[i][j]): 347 self._recursive_split(bbox_partition[i][j], zoom_level + 1, 2 * column + i, 2 * row + 1 - j) 348 349 350class TileSplitter(AreaSplitter): 351 """ A tool that splits the given area into smaller parts. Given the area, time interval and data collection it 352 collects info from Sentinel Hub WFS service about all satellite tiles intersecting the area. For each of them 353 it calculates bounding box and if specified it splits these bounding boxes into smaller bounding boxes. Then 354 it filters out the ones that do not intersect the area. If specified by user it can also reduce the sizes of 355 the remaining bounding boxes to best fit the area. 356 """ 357 def __init__(self, shape_list, crs, time_interval, tile_split_shape=1, data_collection=None, 358 config=None, data_source=None, **kwargs): 359 """ 360 :param shape_list: A list of geometrical shapes describing the area of interest 361 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 362 :param crs: Coordinate reference system of the shapes in `shape_list` 363 :type crs: CRS 364 :param time_interval: Interval with start and end date of the form YYYY-MM-DDThh:mm:ss or YYYY-MM-DD 365 :type time_interval: (str, str) 366 :param tile_split_shape: Parameter that describes the shape in which the satellite tile bounding boxes will be 367 split. It can be a tuple of the form `(n, m)` which means the tile bounding boxes will be 368 split into `n` columns and `m` rows. It can also be a single integer `n` which is the same 369 as `(n, n)`. 370 :type split_shape: int or (int, int) 371 :param data_collection: A satellite data collection 372 :type data_collection: DataCollection 373 :param config: A custom instance of config class to override parameters from the saved configuration. 374 :type config: SHConfig or None 375 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly fit 376 the given area geometry from `shape_list`. 377 :type reduce_bbox_sizes: bool 378 :param data_source: A deprecated alternative of data_collection 379 :type data_source: DataCollection 380 """ 381 super().__init__(shape_list, crs, **kwargs) 382 383 data_collection = DataCollection(handle_deprecated_data_source(data_collection, data_source, 384 default=DataCollection.SENTINEL2_L1C)) 385 if data_collection is DataCollection.DEM: 386 raise ValueError('This splitter does not support splitting area by DEM tiles. Please specify some other ' 387 'DataCollection') 388 389 self.time_interval = time_interval 390 self.tile_split_shape = tile_split_shape 391 self.data_collection = data_collection 392 self.config = config or SHConfig() 393 394 self.tile_dict = None 395 396 self._make_split() 397 398 def _make_split(self): 399 """ This method makes the split 400 """ 401 self.tile_dict = {} 402 403 wfs = WebFeatureService(self.area_bbox, self.time_interval, data_collection=self.data_collection, 404 config=self.config) 405 date_list = wfs.get_dates() 406 geometry_list = wfs.get_geometries() 407 for tile_info, (date, geometry) in zip(wfs, zip(date_list, geometry_list)): 408 tile_name = ''.join(tile_info['properties']['path'].split('/')[4:7]) 409 if tile_name not in self.tile_dict: 410 self.tile_dict[tile_name] = {'bbox': BBox(tile_info['properties']['mbr'], 411 crs=tile_info['properties']['crs']), 412 'times': [], 413 'geometries': []} 414 self.tile_dict[tile_name]['times'].append(date) 415 self.tile_dict[tile_name]['geometries'].append(geometry) 416 417 self.tile_dict = {tile_name: tile_props for tile_name, tile_props in self.tile_dict.items() if 418 self._intersects_area(tile_props['bbox'])} 419 420 self.bbox_list = [] 421 self.info_list = [] 422 423 for tile_name, tile_info in self.tile_dict.items(): 424 tile_bbox = tile_info['bbox'] 425 bbox_splitter = BBoxSplitter([tile_bbox.geometry], tile_bbox.crs, 426 split_shape=self.tile_split_shape) 427 428 for bbox, info in zip(bbox_splitter.get_bbox_list(), bbox_splitter.get_info_list()): 429 if self._intersects_area(bbox): 430 info['tile'] = tile_name 431 432 self.bbox_list.append(bbox) 433 self.info_list.append(info) 434 435 def get_tile_dict(self): 436 """ Returns the dictionary of satellite tiles intersecting the area geometry. For each tile they contain info 437 about their bounding box and lists of acquisitions and geometries 438 439 :return: Dictionary containing info about tiles intersecting the area 440 :rtype: dict 441 """ 442 return self.tile_dict 443 444 445class CustomGridSplitter(AreaSplitter): 446 """ Splitting class which can split according to given custom collection of bounding boxes 447 """ 448 def __init__(self, shape_list, crs, bbox_grid, bbox_split_shape=1, **kwargs): 449 """ 450 :param shape_list: A list of geometrical shapes describing the area of interest 451 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 452 :param crs: Coordinate reference system of the shapes in `shape_list` 453 :type crs: CRS 454 :param bbox_grid: A collection of bounding boxes defining a grid of splitting. All of them have to be in the 455 same CRS. 456 :type bbox_grid: list(BBox) or BBoxCollection 457 :param bbox_split_shape: Parameter that describes the shape in which each of the bounding boxes in the given 458 grid will be split. It can be a tuple of the form `(n, m)` which means the tile bounding boxes will be 459 split into `n` columns and `m` rows. It can also be a single integer `n` which is the same as `(n, n)`. 460 :type bbox_split_shape: int or (int, int) 461 :param reduce_bbox_sizes: If `True` it will reduce the sizes of bounding boxes so that they will tightly fit 462 the given geometry in `shape_list`. 463 :type reduce_bbox_sizes: bool 464 """ 465 super().__init__(shape_list, crs, **kwargs) 466 467 self.bbox_grid = self._parse_bbox_grid(bbox_grid) 468 self.bbox_split_shape = bbox_split_shape 469 470 self._make_split() 471 472 @staticmethod 473 def _parse_bbox_grid(bbox_grid): 474 """ Helper method for parsing bounding box grid. It will try to parse it into `BBoxCollection` 475 """ 476 if isinstance(bbox_grid, BBoxCollection): 477 return bbox_grid 478 479 if isinstance(bbox_grid, list): 480 return BBoxCollection(bbox_grid) 481 482 raise ValueError(f"Parameter 'bbox_grid' should be an instance of {BBoxCollection}") 483 484 def _make_split(self): 485 """ This method makes the split 486 """ 487 self.bbox_list = [] 488 self.info_list = [] 489 490 for grid_idx, grid_bbox in enumerate(self.bbox_grid): 491 if self._intersects_area(grid_bbox): 492 493 bbox_splitter = BBoxSplitter([grid_bbox.geometry], grid_bbox.crs, 494 split_shape=self.bbox_split_shape) 495 496 for bbox, info in zip(bbox_splitter.get_bbox_list(), bbox_splitter.get_info_list()): 497 if self._intersects_area(bbox): 498 info['grid_index'] = grid_idx 499 500 self.bbox_list.append(bbox) 501 self.info_list.append(info) 502 503 504class BaseUtmSplitter(AreaSplitter): 505 """ Base splitter that returns bboxes of fixed size aligned to UTM zones or UTM grid tiles as defined by the MGRS 506 507 The generated bounding box grid will have coordinates in form of 508 `(N * bbox_size_x + offset_x, M * bbox_size_y + offset_y)` 509 """ 510 def __init__(self, shape_list, crs, bbox_size, offset=None): 511 """ 512 :param shape_list: A list of geometrical shapes describing the area of interest 513 :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon) 514 :param crs: Coordinate reference system of the shapes in `shape_list` 515 :type crs: CRS 516 :param bbox_size: A size of generated bounding boxes in horizontal and vertical directions in meters. If a 517 single value is given that will be interpreted as (value, value). 518 :type bbox_size: int or (int, int) or float or (float, float) 519 :param offset: Bounding box offset in horizontal and vertical directions in meters. 520 :type offset: (int, int) or (float, float) or None 521 """ 522 super().__init__(shape_list, crs) 523 524 self.bbox_size = self._parse_split_parameters(bbox_size, allow_float=True) 525 self.offset = self._parse_offset(offset) 526 527 self.shape_geometry = Geometry(self.area_shape, self.crs).transform(CRS.WGS84) 528 529 self.utm_grid = self._get_utm_polygons() 530 531 self._make_split() 532 533 @staticmethod 534 def _parse_offset(offset_input): 535 """ Validates and parses offset input 536 """ 537 if offset_input is None: 538 return 0, 0 539 if isinstance(offset_input, (tuple, list)) and len(offset_input) == 2: 540 return tuple(offset_input) 541 raise ValueError(f'An offset parameter should be a tuple of two numbers, instead {offset_input} was given') 542 543 @abstractmethod 544 def _get_utm_polygons(self): 545 raise NotImplementedError 546 547 @staticmethod 548 def _get_utm_from_props(utm_dict): 549 """ Return the UTM CRS corresponding to the UTM described by the properties dictionary 550 551 :param utm_dict: Dictionary reporting name of the UTM zone and MGRS grid 552 :type utm_dict: dict 553 :return: UTM coordinate reference system 554 :rtype: sentinelhub.CRS 555 """ 556 hemisphere_digit = 6 if utm_dict['direction'] == 'N' else 7 557 zone_number = utm_dict['zone'] 558 return CRS(f'32{hemisphere_digit}{zone_number:02d}') 559 560 def _align_bbox_to_size(self, bbox): 561 """ Align input bbox coordinates to be multiples of the bbox size 562 563 :param bbox: Bounding box in UTM coordinates 564 :type bbox: sentinelhub.BBox 565 :return: BBox objects with coordinates multiples of the bbox size 566 :rtype: sentinelhub.BBox 567 """ 568 size_x, size_y = self.bbox_size 569 offset_x, offset_y = self.offset 570 lower_left_x, lower_left_y = bbox.lower_left 571 572 aligned_x = math.floor((lower_left_x - offset_x) / size_x) * size_x + offset_x 573 aligned_y = math.floor((lower_left_y - offset_y) / size_y) * size_y + offset_y 574 575 return BBox(((aligned_x, aligned_y), bbox.upper_right), crs=bbox.crs) 576 577 def _make_split(self): 578 """ Split each UTM grid into equally sized bboxes in correct UTM zone 579 """ 580 size_x, size_y = self.bbox_size 581 self.bbox_list = [] 582 self.info_list = [] 583 584 index = 0 585 586 for utm_cell in self.utm_grid: 587 utm_cell_geom, utm_cell_prop = utm_cell 588 # the UTM MGRS grid definition contains four 0 zones at the poles (0A, 0B, 0Y, 0Z) 589 if utm_cell_prop['zone'] == 0: 590 continue 591 utm_crs = self._get_utm_from_props(utm_cell_prop) 592 593 intersection = utm_cell_geom.intersection(self.shape_geometry.geometry) 594 595 if not intersection.is_empty and isinstance(intersection, GeometryCollection): 596 intersection = MultiPolygon(geo_object for geo_object in intersection 597 if isinstance(geo_object, (Polygon, MultiPolygon))) 598 599 if not intersection.is_empty: 600 intersection = Geometry(intersection, CRS.WGS84).transform(utm_crs) 601 602 bbox_partition = self._align_bbox_to_size(intersection.bbox).get_partition(size_x=size_x, size_y=size_y) 603 604 columns, rows = len(bbox_partition), len(bbox_partition[0]) 605 for i, j in itertools.product(range(columns), range(rows)): 606 if bbox_partition[i][j].geometry.intersects(intersection.geometry): 607 self.bbox_list.append(bbox_partition[i][j]) 608 self.info_list.append(dict(crs=utm_crs.name, 609 utm_zone=str(utm_cell_prop['zone']).zfill(2), 610 utm_row=utm_cell_prop['row'], 611 direction=utm_cell_prop['direction'], 612 index=index, 613 index_x=i, 614 index_y=j)) 615 index += 1 616 617 def get_bbox_list(self, buffer=None): 618 """ Get list of bounding boxes. 619 620 The CRS is fixed to the computed UTM CRS. This BBox splitter does not support reducing size of output 621 bounding boxes 622 623 :param buffer: A percentage of each BBox size increase. This will cause neighbouring bounding boxes to overlap. 624 :type buffer: float or None 625 :return: List of bounding boxes 626 :rtype: list(BBox) 627 """ 628 return super().get_bbox_list(buffer=buffer) 629 630 631class UtmGridSplitter(BaseUtmSplitter): 632 """ Splitter that returns bounding boxes of fixed size aligned to the UTM MGRS grid 633 """ 634 def _get_utm_polygons(self): 635 """ Find UTM grid zones overlapping with input area shape 636 637 :return: List of geometries and properties of UTM grid zones overlapping with input area shape 638 :rtype: list 639 """ 640 # file downloaded from faculty.baruch.cuny.edu/geoportal/data/esri/world/utmzone.zip 641 utm_grid_filename = os.path.join(os.path.dirname(__file__), '.utmzones.geojson') 642 643 if not os.path.isfile(utm_grid_filename): 644 raise IOError(f'UTM grid definition file does not exist: {os.path.abspath(utm_grid_filename)}') 645 646 with open(utm_grid_filename) as utm_grid_file: 647 utm_grid = json.load(utm_grid_file)['features'] 648 649 utm_geom_list = [shapely.geometry.shape(utm_zone['geometry']) for utm_zone in utm_grid] 650 utm_prop_list = [dict(zone=utm_zone['properties']['ZONE'], 651 row=utm_zone['properties']['ROW_'], 652 direction='N' if utm_zone['properties']['ROW_'] >= 'N' else 'S') for utm_zone in utm_grid] 653 654 return list(zip(utm_geom_list, utm_prop_list)) 655 656 657class UtmZoneSplitter(BaseUtmSplitter): 658 """ Splitter that returns bounding boxes of fixed size aligned to the equator and the UTM zones. 659 """ 660 LNG_MIN, LNG_MAX, LNG_UTM = -180, 180, 6 661 LAT_MIN, LAT_MAX, LAT_EQ = -80, 84, 0 662 663 def _get_utm_polygons(self): 664 """ Find UTM zones overlapping with input area shape 665 666 The returned geometry corresponds to the a triangle ranging from the equator to the north/south pole 667 668 :return: List of geometries and properties of UTM zones overlapping with input area shape 669 :rtype: list 670 """ 671 utm_geom_list = [] 672 for lat in [(self.LAT_EQ, self.LAT_MAX), (self.LAT_MIN, self.LAT_EQ)]: 673 for lng in range(self.LNG_MIN, self.LNG_MAX, self.LNG_UTM): 674 points = [] 675 # A new point is added per each degree - this is inline with geometries used by UtmGridSplitter 676 # In the future the number of points will be calculated according to bbox_size parameter 677 for degree in range(lat[0], lat[1]): 678 points.append((lng, degree)) 679 for degree in range(lng, lng + self.LNG_UTM): 680 points.append((degree, lat[1])) 681 for degree in range(lat[1], lat[0], -1): 682 points.append((lng + self.LNG_UTM, degree)) 683 for degree in range(lng + self.LNG_UTM, lng, -1): 684 points.append((degree, lat[0])) 685 686 utm_geom_list.append(Polygon(points)) 687 688 utm_prop_list = [dict(zone=zone, row='', direction=direction) 689 for direction in ['N', 'S'] for zone in range(1, 61)] 690 691 return list(zip(utm_geom_list, utm_prop_list)) 692 693 694class BatchSplitter(AreaSplitter): 695 """ A splitter that obtains split bounding boxes from Sentinel Hub Batch API 696 """ 697 def __init__(self, *, request_id=None, batch_request=None, config=None): 698 """ 699 :param request_id: An ID of a batch request 700 :type request_id: str or None 701 :param batch_request: A batch request object. It is an alternative to the `request_id` parameter 702 :type batch_request: BatchRequest or None 703 :param config: A configuration object with credentials and information about which service deployment to 704 use. 705 :type config: SHConfig or None 706 """ 707 self.batch_client = SentinelHubBatch(config=config) 708 709 if not (request_id or batch_request): 710 raise ValueError('One of the parameters request_id and batch_request has to be given') 711 if batch_request is None: 712 batch_request = self.batch_client.get_request(request_id) 713 714 self.batch_request = batch_request 715 716 batch_geometry = batch_request.geometry 717 super().__init__([batch_geometry.geometry], batch_geometry.crs) 718 719 self._make_split() 720 721 def _make_split(self): 722 """ This method actually loads bounding boxes from the service and prepares the lists 723 """ 724 tile_info_list = list(self.batch_client.iter_tiles(self.batch_request)) 725 726 tile_geometries = [Geometry.from_geojson(tile_info['geometry']) for tile_info in tile_info_list] 727 original_crs_list = [CRS(tile_info['origin']['crs']['properties']['name']) for tile_info in tile_info_list] 728 729 self.bbox_list = [geometry.transform(crs).bbox for geometry, crs in zip(tile_geometries, original_crs_list)] 730 self.info_list = [ 731 {key: value for key, value in tile_info.items() if key != 'geometry'} for tile_info in tile_info_list 732 ] 733