1""" 2Module implementing geometry classes 3""" 4from abc import ABC, abstractmethod 5from math import ceil 6 7import shapely.geometry 8import shapely.ops 9import shapely.wkt 10 11from .constants import CRS 12from .geo_utils import transform_point 13 14 15class BaseGeometry(ABC): 16 """ Base geometry class 17 """ 18 def __init__(self, crs): 19 """ 20 :param crs: Coordinate reference system of the geometry 21 :type crs: constants.CRS 22 """ 23 self._crs = CRS(crs) 24 25 @property 26 def crs(self): 27 """ Returns the coordinate reference system (CRS) 28 29 :return: Coordinate reference system Enum 30 :rtype: constants.CRS 31 """ 32 return self._crs 33 34 @property 35 @abstractmethod 36 def geometry(self): 37 """ An abstract property - ever subclass must implement geometry property 38 """ 39 raise NotImplementedError 40 41 @property 42 def geojson(self): 43 """ Returns representation in a GeoJSON format. Use json.dump for writing it to file. 44 45 :return: A dictionary in GeoJSON format 46 :rtype: dict 47 """ 48 return self.get_geojson(with_crs=True) 49 50 def get_geojson(self, with_crs=True): 51 """ Returns representation in a GeoJSON format. Use json.dump for writing it to file. 52 53 :param with_crs: A flag indicating if GeoJSON dictionary should contain CRS part 54 :type with_crs: bool 55 :return: A dictionary in GeoJSON format 56 :rtype: dict 57 """ 58 geometry_geojson = shapely.geometry.mapping(self.geometry) 59 60 if with_crs: 61 return { 62 **self._crs_to_geojson(), 63 **geometry_geojson 64 } 65 return geometry_geojson 66 67 def _crs_to_geojson(self): 68 """ Helper method which generates part of GeoJSON format related to CRS 69 """ 70 return { 71 'crs': { 72 'type': 'name', 73 'properties': {'name': f'urn:ogc:def:crs:EPSG::{self.crs.value}'} 74 } 75 } 76 77 @property 78 def wkt(self): 79 """ Transforms geometry object into `Well-known text` format 80 81 :return: string in WKT format 82 :rtype: str 83 """ 84 return self.geometry.wkt 85 86 87class BBox(BaseGeometry): 88 """ Class representing a bounding box in a given CRS. 89 90 Throughout the sentinelhub package this class serves as the canonical representation of a bounding 91 box. It can initialize itself from multiple representations: 92 93 1) ``((min_x,min_y),(max_x,max_y))``, 94 2) ``(min_x,min_y,max_x,max_y)``, 95 3) ``[min_x,min_y,max_x,max_y]``, 96 4) ``[[min_x, min_y],[max_x,max_y]]``, 97 5) ``[(min_x, min_y),(max_x,max_y)]``, 98 6) ``([min_x, min_y],[max_x,max_y])``, 99 7) ``'min_x,min_y,max_x,max_y'``, 100 8) ``{'min_x':min_x, 'max_x':max_x, 'min_y':min_y, 'max_y':max_y}``, 101 9) ``bbox``, where ``bbox`` is an instance of ``BBox``. 102 103 Note that BBox coordinate system depends on ``crs`` parameter: 104 105 - In case of ``constants.CRS.WGS84`` axis x represents longitude and axis y represents latitude 106 - In case of ``constants.CRS.POP_WEB`` axis x represents easting and axis y represents northing 107 - In case of ``constants.CRS.UTM_*`` axis x represents easting and axis y represents northing 108 """ 109 def __init__(self, bbox, crs): 110 """ 111 :param bbox: A bbox in any valid representation 112 :param crs: Coordinate reference system of the bounding box 113 :type crs: constants.CRS 114 """ 115 x_fst, y_fst, x_snd, y_snd = BBox._to_tuple(bbox) 116 self.min_x = min(x_fst, x_snd) 117 self.max_x = max(x_fst, x_snd) 118 self.min_y = min(y_fst, y_snd) 119 self.max_y = max(y_fst, y_snd) 120 121 super().__init__(crs) 122 123 def __iter__(self): 124 """ This method enables iteration over coordinates of bounding box 125 """ 126 return iter(self.lower_left + self.upper_right) 127 128 def __repr__(self): 129 """ Class representation 130 """ 131 return f'{self.__class__.__name__}(({self.lower_left}, {self.upper_right}), crs={repr(self.crs)})' 132 133 def __str__(self, reverse=False): 134 """ Transforms bounding box into a string of coordinates 135 136 :param reverse: `True` if x and y coordinates should be switched and `False` otherwise 137 :type reverse: bool 138 :return: String of coordinates 139 :rtype: str 140 """ 141 if reverse: 142 return f'{self.min_y},{self.min_x},{self.max_y},{self.max_x}' 143 return f'{self.min_x},{self.min_y},{self.max_x},{self.max_y}' 144 145 def __eq__(self, other): 146 """ Method for comparing two bounding boxes 147 148 :param other: Another bounding box object 149 :type other: BBox 150 :return: `True` if bounding boxes have the same coordinates and the same CRS and `False otherwise 151 :rtype: bool 152 """ 153 if not isinstance(other, BBox): 154 return False 155 return list(self) == list(other) and self.crs is other.crs 156 157 @property 158 def lower_left(self): 159 """ Returns the lower left vertex of the bounding box 160 161 :return: min_x, min_y 162 :rtype: (float, float) 163 """ 164 return self.min_x, self.min_y 165 166 @property 167 def upper_right(self): 168 """ Returns the upper right vertex of the bounding box 169 170 :return: max_x, max_y 171 :rtype: (float, float) 172 """ 173 return self.max_x, self.max_y 174 175 @property 176 def middle(self): 177 """ Returns the middle point of the bounding box 178 179 :return: middle point 180 :rtype: (float, float) 181 """ 182 return (self.min_x + self.max_x) / 2, (self.min_y + self.max_y) / 2 183 184 def reverse(self): 185 """ Returns a new BBox object where x and y coordinates are switched 186 187 :return: New BBox object with switched coordinates 188 :rtype: BBox 189 """ 190 return BBox((self.min_y, self.min_x, self.max_y, self.max_x), crs=self.crs) 191 192 def transform(self, crs, always_xy=True): 193 """ Transforms BBox from current CRS to target CRS 194 195 This transformation will take lower left and upper right corners of the bounding box, transform these 2 points 196 and define a new bounding box with them. The resulting bounding box might not completely cover the original 197 bounding box but at least the transformation is reversible. 198 199 :param crs: target CRS 200 :type crs: constants.CRS 201 :param always_xy: Parameter that is passed to `pyproj.Transformer` object and defines axis order for 202 transformation. The default value `True` is in most cases the correct one. 203 :type always_xy: bool 204 :return: Bounding box in target CRS 205 :rtype: BBox 206 """ 207 new_crs = CRS(crs) 208 return BBox((transform_point(self.lower_left, self.crs, new_crs, always_xy=always_xy), 209 transform_point(self.upper_right, self.crs, new_crs, always_xy=always_xy)), crs=new_crs) 210 211 def transform_bounds(self, crs, always_xy=True): 212 """ Alternative way to transform BBox from current CRS to target CRS. 213 214 This transformation will transform the bounding box geometry to another CRS as a geometric object, and then 215 define a new bounding box from boundaries of that geometry. The resulting bounding box might be larger than 216 original bounding box but it will always completely cover it. 217 218 :param crs: target CRS 219 :type crs: constants.CRS 220 :param always_xy: Parameter that is passed to `pyproj.Transformer` object and defines axis order for 221 transformation. The default value `True` is in most cases the correct one. 222 :type always_xy: bool 223 :return: Bounding box in target CRS 224 :rtype: BBox 225 """ 226 bbox_geometry = Geometry(self.geometry, self.crs) 227 bbox_geometry = bbox_geometry.transform(crs, always_xy=always_xy) 228 return bbox_geometry.bbox 229 230 def buffer(self, buffer): 231 """ Changes both BBox dimensions (width and height) by a percentage of size of each dimension. If number is 232 negative, the size will decrease. Returns a new instance of BBox object. 233 234 :param buffer: A percentage of BBox size change 235 :type buffer: float 236 :return: A new bounding box of buffered size 237 :rtype: BBox 238 """ 239 if buffer < -1: 240 raise ValueError('Cannot reduce the bounding box to nothing, buffer must be >= -1.0') 241 ratio = 1 + buffer 242 mid_x, mid_y = self.middle 243 return BBox((mid_x - (mid_x - self.min_x) * ratio, mid_y - (mid_y - self.min_y) * ratio, 244 mid_x + (self.max_x - mid_x) * ratio, mid_y + (self.max_y - mid_y) * ratio), self.crs) 245 246 def get_polygon(self, reverse=False): 247 """ Returns a tuple of coordinates of 5 points describing a polygon. Points are listed in clockwise order, first 248 point is the same as the last. 249 250 :param reverse: `True` if x and y coordinates should be switched and `False` otherwise 251 :type reverse: bool 252 :return: `((x_1, y_1), ... , (x_5, y_5))` 253 :rtype: tuple(tuple(float)) 254 """ 255 bbox = self.reverse() if reverse else self 256 polygon = ((bbox.min_x, bbox.min_y), 257 (bbox.min_x, bbox.max_y), 258 (bbox.max_x, bbox.max_y), 259 (bbox.max_x, bbox.min_y), 260 (bbox.min_x, bbox.min_y)) 261 return polygon 262 263 @property 264 def geometry(self): 265 """ Returns polygon geometry in shapely format 266 267 :return: A polygon in shapely format 268 :rtype: shapely.geometry.polygon.Polygon 269 """ 270 return shapely.geometry.Polygon(self.get_polygon()) 271 272 def get_partition(self, num_x=None, num_y=None, size_x=None, size_y=None): 273 """ Partitions bounding box into smaller bounding boxes of the same size. 274 275 If `num_x` and `num_y` are specified, the total number of BBoxes is know but not the size. If `size_x` and 276 `size_y` are provided, the BBox size is fixed but the number of BBoxes is not known in advance. In the latter 277 case, the generated bounding boxes might cover an area larger than the parent BBox. 278 279 :param num_x: Number of parts BBox will be horizontally divided into. 280 :type num_x: int or None 281 :param num_y: Number of parts BBox will be vertically divided into. 282 :type num_y: int or None 283 :param size_x: Physical dimension of BBox along easting coordinate 284 :type size_x: float or None 285 :param size_y: Physical dimension of BBox along northing coordinate 286 :type size_y: float or None 287 :return: Two-dimensional list of smaller bounding boxes. Their location is 288 :rtype: list(list(BBox)) 289 """ 290 if (num_x is not None and num_y is not None) and (size_x is None and size_y is None): 291 size_x, size_y = (self.max_x - self.min_x) / num_x, (self.max_y - self.min_y) / num_y 292 elif (size_x is not None and size_y is not None) and (num_x is None and num_y is None): 293 num_x, num_y = ceil((self.max_x - self.min_x) / size_x), ceil((self.max_y - self.min_y) / size_y) 294 else: 295 raise ValueError('Not supported partition. Either (num_x, num_y) or (size_x, size_y) must be specified') 296 297 return [[BBox([self.min_x + i * size_x, self.min_y + j * size_y, 298 self.min_x + (i + 1) * size_x, self.min_y + (j + 1) * size_y], 299 crs=self.crs) for j in range(num_y)] for i in range(num_x)] 300 301 def get_transform_vector(self, resx, resy): 302 """ Given resolution it returns a transformation vector 303 304 :param resx: Resolution in x direction 305 :type resx: float or int 306 :param resy: Resolution in y direction 307 :type resy: float or int 308 :return: A tuple with 6 numbers representing transformation vector 309 :rtype: tuple(float) 310 """ 311 return self.min_x, self._parse_resolution(resx), 0, self.max_y, 0, -self._parse_resolution(resy) 312 313 @staticmethod 314 def _parse_resolution(res): 315 """ Helper method for parsing given resolution. It will also try to parse a string into float 316 317 :return: A float value of resolution 318 :rtype: float 319 """ 320 if isinstance(res, str): 321 return float(res.strip('m')) 322 if isinstance(res, (int, float)): 323 return float(res) 324 325 raise TypeError(f'Resolution should be a float, got resolution of type {type(res)}') 326 327 @staticmethod 328 def _to_tuple(bbox): 329 """ Converts the input bbox representation (see the constructor docstring for a list of valid representations) 330 into a flat tuple 331 332 :param bbox: A bbox in one of 7 forms listed in the class description. 333 :return: A flat tuple of size 334 :raises: TypeError 335 """ 336 if isinstance(bbox, (list, tuple)): 337 return BBox._tuple_from_list_or_tuple(bbox) 338 if isinstance(bbox, str): 339 return BBox._tuple_from_str(bbox) 340 if isinstance(bbox, dict): 341 return BBox._tuple_from_dict(bbox) 342 if isinstance(bbox, BBox): 343 return BBox._tuple_from_bbox(bbox) 344 if isinstance(bbox, shapely.geometry.base.BaseGeometry): 345 return bbox.bounds 346 raise TypeError('Invalid bbox representation') 347 348 @staticmethod 349 def _tuple_from_list_or_tuple(bbox): 350 """ Converts a list or tuple representation of a bbox into a flat tuple representation. 351 352 :param bbox: a list or tuple with 4 coordinates that is either flat or nested 353 :return: tuple (min_x,min_y,max_x,max_y) 354 :raises: TypeError 355 """ 356 if len(bbox) == 4: 357 return tuple(map(float, bbox)) 358 if len(bbox) == 2 and all(isinstance(point, (list, tuple)) for point in bbox): 359 return BBox._tuple_from_list_or_tuple(bbox[0] + bbox[1]) 360 raise TypeError('Expected a valid list or tuple representation of a bbox') 361 362 @staticmethod 363 def _tuple_from_str(bbox): 364 """ Parses a string of numbers separated by any combination of commas and spaces 365 366 :param bbox: e.g. str of the form `min_x ,min_y max_x, max_y` 367 :return: tuple (min_x,min_y,max_x,max_y) 368 """ 369 return tuple(float(s) for s in bbox.replace(',', ' ').split() if s) 370 371 @staticmethod 372 def _tuple_from_dict(bbox): 373 """ Converts a dictionary representation of a bbox into a flat tuple representation 374 375 :param bbox: a dict with keys "min_x, "min_y", "max_x", and "max_y" 376 :return: tuple (min_x,min_y,max_x,max_y) 377 :raises: KeyError 378 """ 379 return bbox['min_x'], bbox['min_y'], bbox['max_x'], bbox['max_y'] 380 381 @staticmethod 382 def _tuple_from_bbox(bbox): 383 """ Converts a BBox instance into a tuple 384 385 :param bbox: An instance of the BBox type 386 :return: tuple (min_x, min_y, max_x, max_y) 387 """ 388 return bbox.lower_left + bbox.upper_right 389 390 391class Geometry(BaseGeometry): 392 """ A class that combines shapely geometry with coordinate reference system. It currently supports polygons and 393 multipolygons. 394 395 It can be initialize with any of the following geometry representations: 396 - `shapely.geometry.Polygon` or `shapely.geometry.MultiPolygon` 397 - A GeoJSON dictionary with (multi)polygon coordinates 398 - A WKT string with (multi)polygon coordinates 399 """ 400 def __init__(self, geometry, crs): 401 """ 402 :param geometry: A polygon or multipolygon in any valid representation 403 :type geometry: shapely.geometry.Polygon or shapely.geometry.MultiPolygon or dict or str 404 :param crs: Coordinate reference system of the geometry 405 :type crs: constants.CRS 406 """ 407 self._geometry = self._parse_geometry(geometry) 408 409 super().__init__(crs) 410 411 def __repr__(self): 412 """ Method for class representation 413 """ 414 return f'{self.__class__.__name__}({self.wkt}, crs={repr(self.crs)})' 415 416 def __eq__(self, other): 417 """ Method for comparing two Geometry classes 418 419 :param other: Another Geometry object 420 :type other: Geometry 421 :return: `True` if geometry objects have the same geometry and CRS and `False` otherwise 422 :rtype: bool 423 """ 424 if not isinstance(other, Geometry): 425 return False 426 return self.geometry == other.geometry and self.crs is other.crs 427 428 def reverse(self): 429 """ Returns a new Geometry object where x and y coordinates are switched 430 431 :return: New Geometry object with switched coordinates 432 :rtype: Geometry 433 """ 434 return Geometry(shapely.ops.transform(lambda x, y: (y, x), self.geometry), crs=self.crs) 435 436 def transform(self, crs, always_xy=True): 437 """ Transforms Geometry from current CRS to target CRS 438 439 :param crs: target CRS 440 :type crs: constants.CRS 441 :param always_xy: Parameter that is passed to `pyproj.Transformer` object and defines axis order for 442 transformation. The default value `True` is in most cases the correct one. 443 :type always_xy: bool 444 :return: Geometry in target CRS 445 :rtype: Geometry 446 """ 447 new_crs = CRS(crs) 448 449 geometry = self.geometry 450 if new_crs is not self.crs: 451 transform_function = self.crs.get_transform_function(new_crs, always_xy=always_xy) 452 geometry = shapely.ops.transform(transform_function, geometry) 453 454 return Geometry(geometry, crs=new_crs) 455 456 @classmethod 457 def from_geojson(cls, geojson, crs=None): 458 """ Create Geometry object from geojson. It will parse crs from geojson (if info is available), 459 otherwise it will be set to crs (WGS84 if parameter is empty) 460 461 :param geojson: geojson geometry (single feature) 462 :param crs: crs to be used if not available in geojson, CRS.WGS84 if not provided 463 :return: Geometry object 464 """ 465 try: 466 crs = CRS(geojson['crs']['properties']['name']) 467 except (KeyError, AttributeError, TypeError): 468 pass 469 470 if not crs: 471 crs = CRS.WGS84 472 473 return cls(geojson, crs=crs) 474 475 @property 476 def geometry(self): 477 """ Returns shapely object representing geometry in this class 478 479 :return: A polygon or a multipolygon in shapely format 480 :rtype: shapely.geometry.Polygon or shapely.geometry.MultiPolygon 481 """ 482 return self._geometry 483 484 @property 485 def bbox(self): 486 """ Returns BBox object representing bounding box around the geometry 487 488 :return: A bounding box, with same CRS 489 :rtype: BBox 490 """ 491 return BBox(self.geometry, self.crs) 492 493 @staticmethod 494 def _parse_geometry(geometry): 495 """ Parses given geometry into shapely object 496 497 :param geometry: 498 :return: Shapely polygon or multipolygon 499 :rtype: shapely.geometry.Polygon or shapely.geometry.MultiPolygon 500 :raises TypeError 501 """ 502 if isinstance(geometry, str): 503 geometry = shapely.wkt.loads(geometry) 504 elif isinstance(geometry, dict): 505 geometry = shapely.geometry.shape(geometry) 506 elif not isinstance(geometry, shapely.geometry.base.BaseGeometry): 507 raise TypeError('Unsupported geometry representation') 508 509 if not isinstance(geometry, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon)): 510 raise ValueError(f'Supported geometry types are polygon and multipolygon, got {type(geometry)}') 511 512 return geometry 513 514 515class BBoxCollection(BaseGeometry): 516 """ A collection of bounding boxes 517 """ 518 def __init__(self, bbox_list): 519 """ 520 :param bbox_list: A list of BBox objects which have to be in the same CRS 521 :type bbox_list: list(BBox) 522 """ 523 self._bbox_list, crs = self._parse_bbox_list(bbox_list) 524 self._geometry = self._get_geometry() 525 526 super().__init__(crs) 527 528 def __repr__(self): 529 """ Method for class representation 530 """ 531 bbox_list_repr = ', '.join([repr(bbox) for bbox in self.bbox_list]) 532 return f'{self.__class__.__name__}({bbox_list_repr})' 533 534 def __eq__(self, other): 535 """ Method for comparing two BBoxCollection classes 536 """ 537 if not isinstance(other, BBoxCollection): 538 return False 539 return self.crs is other.crs and len(self.bbox_list) == len(other.bbox_list) and \ 540 all(bbox == other_bbox for bbox, other_bbox in zip(self, other)) 541 542 def __iter__(self): 543 """ This method enables iteration over bounding boxes in collection 544 """ 545 return iter(self.bbox_list) 546 547 @property 548 def bbox_list(self): 549 """ Returns the list of bounding boxes from collection 550 551 :return: The list of bounding boxes 552 :rtype: list(BBox) 553 """ 554 return self._bbox_list 555 556 @property 557 def geometry(self): 558 """ Returns shapely object representing geometry 559 560 :return: A multipolygon of bounding boxes 561 :rtype: shapely.geometry.MultiPolygon 562 """ 563 return self._geometry 564 565 @property 566 def bbox(self): 567 """ Returns BBox object representing bounding box around the geometry 568 569 :return: A bounding box, with same CRS 570 :rtype: BBox 571 """ 572 return BBox(self.geometry, self.crs) 573 574 def reverse(self): 575 """ Returns a new BBoxCollection object where all x and y coordinates are switched 576 577 :return: New Geometry object with switched coordinates 578 :rtype: BBoxCollection 579 """ 580 return BBoxCollection([bbox.reverse() for bbox in self.bbox_list]) 581 582 def transform(self, crs, always_xy=True): 583 """ Transforms BBoxCollection from current CRS to target CRS 584 585 :param crs: target CRS 586 :type crs: constants.CRS 587 :param always_xy: Parameter that is passed to `pyproj.Transformer` object and defines axis order for 588 transformation. The default value `True` is in most cases the correct one. 589 :type always_xy: bool 590 :return: BBoxCollection in target CRS 591 :rtype: BBoxCollection 592 """ 593 return BBoxCollection([bbox.transform(crs, always_xy=always_xy) for bbox in self.bbox_list]) 594 595 def _get_geometry(self): 596 """ Creates a multipolygon of bounding box polygons 597 """ 598 return shapely.geometry.MultiPolygon([bbox.geometry for bbox in self.bbox_list]) 599 600 @staticmethod 601 def _parse_bbox_list(bbox_list): 602 """ Helper method for parsing a list of bounding boxes 603 """ 604 if isinstance(bbox_list, BBoxCollection): 605 return bbox_list.bbox_list, bbox_list.crs 606 607 if not isinstance(bbox_list, list) or not bbox_list: 608 raise ValueError('Expected non-empty list of BBox objects') 609 610 for bbox in bbox_list: 611 if not isinstance(bbox, BBox): 612 raise ValueError(f'Elements in the list should be of type {BBox.__name__}, got {type(bbox)}') 613 614 crs = bbox_list[0].crs 615 for bbox in bbox_list: 616 if bbox.crs is not crs: 617 raise ValueError('All bounding boxes should have the same CRS') 618 619 return bbox_list, crs 620