1"""Base geometry class and utilities 2 3Note: a third, z, coordinate value may be used when constructing 4geometry objects, but has no effect on geometric analysis. All 5operations are performed in the x-y plane. Thus, geometries with 6different z values may intersect or be equal. 7""" 8 9from binascii import a2b_hex 10from ctypes import pointer, c_size_t, c_char_p, c_void_p 11from itertools import islice 12import logging 13import math 14import sys 15from warnings import warn 16from functools import wraps 17import warnings 18 19from shapely.affinity import affine_transform 20from shapely.coords import CoordinateSequence 21from shapely.errors import GeometryTypeError, WKBReadingError, WKTReadingError 22from shapely.errors import ShapelyDeprecationWarning 23from shapely.geos import WKBWriter, WKTWriter 24from shapely.geos import lgeos 25from shapely.impl import DefaultImplementation, delegated 26 27log = logging.getLogger(__name__) 28 29try: 30 import numpy as np 31 integer_types = (int, np.integer) 32except ImportError: 33 integer_types = (int,) 34 35 36GEOMETRY_TYPES = [ 37 'Point', 38 'LineString', 39 'LinearRing', 40 'Polygon', 41 'MultiPoint', 42 'MultiLineString', 43 'MultiPolygon', 44 'GeometryCollection', 45] 46 47 48def dump_coords(geom): 49 """Dump coordinates of a geometry in the same order as data packing""" 50 if not isinstance(geom, BaseGeometry): 51 raise ValueError('Must be instance of a geometry class; found ' + 52 geom.__class__.__name__) 53 elif geom.type in ('Point', 'LineString', 'LinearRing'): 54 return geom.coords[:] 55 elif geom.type == 'Polygon': 56 return geom.exterior.coords[:] + [i.coords[:] for i in geom.interiors] 57 elif geom.type.startswith('Multi') or geom.type == 'GeometryCollection': 58 # Recursive call 59 return [dump_coords(part) for part in geom.geoms] 60 else: 61 raise GeometryTypeError('Unhandled geometry type: ' + repr(geom.type)) 62 63 64def geometry_type_name(g): 65 if g is None: 66 raise ValueError("Null geometry has no type") 67 return GEOMETRY_TYPES[lgeos.GEOSGeomTypeId(g)] 68 69 70def geom_factory(g, parent=None): 71 # Abstract geometry factory for use with topological methods below 72 if not g: 73 raise ValueError("No Shapely geometry can be created from null value") 74 ob = BaseGeometry() 75 geom_type = geometry_type_name(g) 76 # TODO: check cost of dynamic import by profiling 77 mod = __import__( 78 'shapely.geometry', 79 globals(), 80 locals(), 81 [geom_type], 82 ) 83 ob.__class__ = getattr(mod, geom_type) 84 ob._set_geom(g) 85 ob.__p__ = parent 86 if lgeos.methods['has_z'](g): 87 ob._ndim = 3 88 else: 89 ob._ndim = 2 90 ob._is_empty = False 91 return ob 92 93 94def deserialize_wkb(data): 95 geom = lgeos.GEOSGeomFromWKB_buf(c_char_p(data), c_size_t(len(data))) 96 if not geom: 97 raise WKBReadingError( 98 "Could not create geometry because of errors while reading input.") 99 return geom 100 101 102def geos_geom_from_py(ob, create_func=None): 103 """Helper function for geos_*_from_py functions in each geom type. 104 105 If a create_func is specified the coordinate sequence is cloned and a new 106 geometry is created with it, otherwise the geometry is cloned directly. 107 This behaviour is useful for converting between LineString and LinearRing 108 objects. 109 """ 110 if create_func is None: 111 geom = lgeos.GEOSGeom_clone(ob._geom) 112 else: 113 cs = lgeos.GEOSGeom_getCoordSeq(ob._geom) 114 cs = lgeos.GEOSCoordSeq_clone(cs) 115 geom = create_func(cs) 116 117 N = ob._ndim 118 119 return geom, N 120 121 122def exceptNull(func): 123 """Decorator which helps avoid GEOS operations on null pointers.""" 124 @wraps(func) 125 def wrapper(*args, **kwargs): 126 if not args[0]._geom or args[0].is_empty: 127 raise ValueError("Null/empty geometry supports no operations") 128 return func(*args, **kwargs) 129 return wrapper 130 131 132class CAP_STYLE: 133 round = 1 134 flat = 2 135 square = 3 136 137 138class JOIN_STYLE: 139 round = 1 140 mitre = 2 141 bevel = 3 142 143EMPTY = deserialize_wkb(a2b_hex(b'010700000000000000')) 144 145 146class BaseGeometry: 147 """ 148 Provides GEOS spatial predicates and topological operations. 149 150 """ 151 152 # Attributes 153 # ---------- 154 # __geom__ : c_void_p 155 # Cached ctypes pointer to GEOS geometry. Not to be accessed. 156 # _geom : c_void_p 157 # Property by which the GEOS geometry is accessed. 158 # __p__ : object 159 # Parent (Shapely) geometry 160 # _ctypes_data : object 161 # Cached ctypes data buffer 162 # _ndim : int 163 # Number of dimensions (2 or 3, generally) 164 # _crs : object 165 # Coordinate reference system. Available for Shapely extensions, but 166 # not implemented here. 167 # _other_owned : bool 168 # True if this object's GEOS geometry is owned by another as in the 169 # case of a multipart geometry member. 170 __geom__ = EMPTY 171 __p__ = None 172 _ctypes_data = None 173 _ndim = None 174 _crs = None 175 _other_owned = False 176 _is_empty = True 177 178 # Backend config 179 impl = DefaultImplementation 180 181 # a reference to the so/dll proxy to preserve access during clean up 182 _lgeos = lgeos 183 184 def empty(self, val=EMPTY): 185 warn( 186 "The 'empty()' method is deprecated and will be removed in " 187 "Shapely 2.0", 188 ShapelyDeprecationWarning, stacklevel=2) 189 self._empty(val=val) 190 191 def _empty(self, val=EMPTY): 192 if not self._other_owned and self.__geom__ and self.__geom__ != EMPTY: 193 try: 194 self._lgeos.GEOSGeom_destroy(self.__geom__) 195 except (AttributeError, TypeError): 196 # _lgeos might be empty on shutdown 197 log.exception("Failed to delete GEOS geom") 198 199 self._is_empty = True 200 self.__geom__ = val 201 202 def __bool__(self): 203 return self.is_empty is False 204 205 def __nonzero__(self): 206 return self.__bool__() 207 208 def __del__(self): 209 self._empty(val=None) 210 self.__p__ = None 211 212 def __str__(self): 213 return self.wkt 214 215 # To support pickling 216 def __reduce__(self): 217 return (self.__class__, (), self.wkb) 218 219 def __setstate__(self, state): 220 self._empty() 221 self.__geom__ = deserialize_wkb(state) 222 self._is_empty = False 223 if lgeos.methods['has_z'](self.__geom__): 224 self._ndim = 3 225 else: 226 self._ndim = 2 227 228 @property 229 def _geom(self): 230 return self.__geom__ 231 232 @_geom.setter 233 def _geom(self, val): 234 warn( 235 "Setting the '_geom' attribute directly is deprecated, and will " 236 "no longer be possible in Shapely 2.0.", 237 ShapelyDeprecationWarning, stacklevel=2) 238 self._set_geom(val) 239 240 def _set_geom(self, val): 241 self._empty() 242 self._is_empty = val in [EMPTY, None] 243 self.__geom__ = val 244 245 def __setattr__(self, name, value): 246 # first try regular attribute access via __getattribute__, so that 247 # our own (existing) attributes don't raise a warning 248 try: 249 object.__getattribute__(self, name) 250 super().__setattr__(name, value) 251 return 252 except AttributeError: 253 pass 254 255 # if custom atribute, raise the warning 256 warn( 257 "Setting custom attributes on geometry objects is deprecated, " 258 "and will raise an AttributeError in Shapely 2.0", 259 ShapelyDeprecationWarning, stacklevel=2 260 ) 261 super().__setattr__(name, value) 262 263 # Operators 264 # --------- 265 266 def __and__(self, other): 267 return self.intersection(other) 268 269 def __or__(self, other): 270 return self.union(other) 271 272 def __sub__(self, other): 273 return self.difference(other) 274 275 def __xor__(self, other): 276 return self.symmetric_difference(other) 277 278 def __eq__(self, other): 279 return ( 280 type(other) == type(self) and 281 tuple(self.coords) == tuple(other.coords) 282 ) 283 284 def __ne__(self, other): 285 return not self.__eq__(other) 286 287 __hash__ = None 288 289 # Array and ctypes interfaces 290 # --------------------------- 291 292 @property 293 def _ctypes(self): 294 raise NotImplementedError 295 296 @property 297 def ctypes(self): 298 """Return ctypes buffer""" 299 warn( 300 "Accessing the 'ctypes' attribute is deprecated," 301 " and will not be possible any more in Shapely 2.0", 302 ShapelyDeprecationWarning, stacklevel=2) 303 return self._ctypes 304 305 @property 306 def _array_interface_base(self): 307 if sys.byteorder == 'little': 308 typestr = '<f8' 309 elif sys.byteorder == 'big': 310 typestr = '>f8' 311 else: 312 raise ValueError( 313 "Unsupported byteorder: neither little nor big-endian") 314 return { 315 'version': 3, 316 'typestr': typestr, 317 'data': self._ctypes, 318 } 319 320 @property 321 def array_interface_base(self): 322 warn( 323 "The 'array_interface_base' property is deprecated and will be " 324 "removed in Shapely 2.0.", 325 ShapelyDeprecationWarning, stacklevel=2) 326 return self._array_interface_base() 327 328 @property 329 def __array_interface__(self): 330 """Provide the Numpy array protocol.""" 331 raise AttributeError("Not implemented") 332 333 # Coordinate access 334 # ----------------- 335 336 def _get_coords(self): 337 """Access to geometry's coordinates (CoordinateSequence)""" 338 raise NotImplementedError("set_coords must be provided by derived classes") 339 340 def _set_coords(self, ob): 341 raise NotImplementedError( 342 "set_coords must be provided by derived classes") 343 344 coords = property(_get_coords, _set_coords) 345 346 @property 347 def xy(self): 348 """Separate arrays of X and Y coordinate values""" 349 raise NotImplementedError 350 351 # Python feature protocol 352 353 @property 354 def __geo_interface__(self): 355 """Dictionary representation of the geometry""" 356 raise NotImplementedError 357 358 # Type of geometry and its representations 359 # ---------------------------------------- 360 361 def geometryType(self): 362 return geometry_type_name(self._geom) 363 364 @property 365 def type(self): 366 return self.geometryType() 367 368 @property 369 def wkt(self): 370 """WKT representation of the geometry""" 371 return WKTWriter(lgeos).write(self) 372 373 @property 374 def wkb(self): 375 """WKB representation of the geometry""" 376 return WKBWriter(lgeos).write(self) 377 378 @property 379 def wkb_hex(self): 380 """WKB hex representation of the geometry""" 381 return WKBWriter(lgeos).write_hex(self) 382 383 def svg(self, scale_factor=1., **kwargs): 384 """Raises NotImplementedError""" 385 raise NotImplementedError 386 387 def _repr_svg_(self): 388 """SVG representation for iPython notebook""" 389 svg_top = '<svg xmlns="http://www.w3.org/2000/svg" ' \ 390 'xmlns:xlink="http://www.w3.org/1999/xlink" ' 391 if self.is_empty: 392 return svg_top + '/>' 393 else: 394 # Establish SVG canvas that will fit all the data + small space 395 xmin, ymin, xmax, ymax = self.bounds 396 if xmin == xmax and ymin == ymax: 397 # This is a point; buffer using an arbitrary size 398 xmin, ymin, xmax, ymax = self.buffer(1).bounds 399 else: 400 # Expand bounds by a fraction of the data ranges 401 expand = 0.04 # or 4%, same as R plots 402 widest_part = max([xmax - xmin, ymax - ymin]) 403 expand_amount = widest_part * expand 404 xmin -= expand_amount 405 ymin -= expand_amount 406 xmax += expand_amount 407 ymax += expand_amount 408 dx = xmax - xmin 409 dy = ymax - ymin 410 width = min([max([100., dx]), 300]) 411 height = min([max([100., dy]), 300]) 412 try: 413 scale_factor = max([dx, dy]) / max([width, height]) 414 except ZeroDivisionError: 415 scale_factor = 1. 416 view_box = "{} {} {} {}".format(xmin, ymin, dx, dy) 417 transform = "matrix(1,0,0,-1,0,{})".format(ymax + ymin) 418 return svg_top + ( 419 'width="{1}" height="{2}" viewBox="{0}" ' 420 'preserveAspectRatio="xMinYMin meet">' 421 '<g transform="{3}">{4}</g></svg>' 422 ).format(view_box, width, height, transform, 423 self.svg(scale_factor)) 424 425 @property 426 def geom_type(self): 427 """Name of the geometry's type, such as 'Point'""" 428 return self.geometryType() 429 430 # Real-valued properties and methods 431 # ---------------------------------- 432 433 @property 434 def area(self): 435 """Unitless area of the geometry (float)""" 436 return self.impl['area'](self) 437 438 def distance(self, other): 439 """Unitless distance to other geometry (float)""" 440 return self.impl['distance'](self, other) 441 442 def hausdorff_distance(self, other): 443 """Unitless hausdorff distance to other geometry (float)""" 444 return self.impl['hausdorff_distance'](self, other) 445 446 @property 447 def length(self): 448 """Unitless length of the geometry (float)""" 449 return self.impl['length'](self) 450 451 @property 452 def minimum_clearance(self): 453 """Unitless distance by which a node could be moved to produce an invalid geometry (float)""" 454 return self.impl['minimum_clearance'](self) 455 456 # Topological properties 457 # ---------------------- 458 459 @property 460 def boundary(self): 461 """Returns a lower dimension geometry that bounds the object 462 463 The boundary of a polygon is a line, the boundary of a line is a 464 collection of points. The boundary of a point is an empty (null) 465 collection. 466 """ 467 return geom_factory(self.impl['boundary'](self)) 468 469 @property 470 def bounds(self): 471 """Returns minimum bounding region (minx, miny, maxx, maxy)""" 472 if self.is_empty: 473 return () 474 else: 475 return self.impl['bounds'](self) 476 477 @property 478 def centroid(self): 479 """Returns the geometric center of the object""" 480 return geom_factory(self.impl['centroid'](self)) 481 482 @delegated 483 def representative_point(self): 484 """Returns a point guaranteed to be within the object, cheaply.""" 485 return geom_factory(self.impl['representative_point'](self)) 486 487 @property 488 def convex_hull(self): 489 """Imagine an elastic band stretched around the geometry: that's a 490 convex hull, more or less 491 492 The convex hull of a three member multipoint, for example, is a 493 triangular polygon. 494 """ 495 return geom_factory(self.impl['convex_hull'](self)) 496 497 @property 498 def envelope(self): 499 """A figure that envelopes the geometry""" 500 return geom_factory(self.impl['envelope'](self)) 501 502 @property 503 def minimum_rotated_rectangle(self): 504 """Returns the general minimum bounding rectangle of 505 the geometry. Can possibly be rotated. If the convex hull 506 of the object is a degenerate (line or point) this same degenerate 507 is returned. 508 """ 509 # first compute the convex hull 510 hull = self.convex_hull 511 try: 512 coords = hull.exterior.coords 513 except AttributeError: # may be a Point or a LineString 514 return hull 515 # generate the edge vectors between the convex hull's coords 516 edges = ((pt2[0] - pt1[0], pt2[1] - pt1[1]) for pt1, pt2 in zip( 517 coords, islice(coords, 1, None))) 518 519 def _transformed_rects(): 520 for dx, dy in edges: 521 # compute the normalized direction vector of the edge 522 # vector. 523 length = math.sqrt(dx ** 2 + dy ** 2) 524 ux, uy = dx / length, dy / length 525 # compute the normalized perpendicular vector 526 vx, vy = -uy, ux 527 # transform hull from the original coordinate system to 528 # the coordinate system defined by the edge and compute 529 # the axes-parallel bounding rectangle. 530 transf_rect = affine_transform( 531 hull, (ux, uy, vx, vy, 0, 0)).envelope 532 # yield the transformed rectangle and a matrix to 533 # transform it back to the original coordinate system. 534 yield (transf_rect, (ux, vx, uy, vy, 0, 0)) 535 536 # check for the minimum area rectangle and return it 537 transf_rect, inv_matrix = min( 538 _transformed_rects(), key=lambda r: r[0].area) 539 return affine_transform(transf_rect, inv_matrix) 540 541 def buffer(self, distance, resolution=16, quadsegs=None, 542 cap_style=CAP_STYLE.round, join_style=JOIN_STYLE.round, 543 mitre_limit=5.0, single_sided=False): 544 """Get a geometry that represents all points within a distance 545 of this geometry. 546 547 A positive distance produces a dilation, a negative distance an 548 erosion. A very small or zero distance may sometimes be used to 549 "tidy" a polygon. 550 551 Parameters 552 ---------- 553 distance : float 554 The distance to buffer around the object. 555 resolution : int, optional 556 The resolution of the buffer around each vertex of the 557 object. 558 quadsegs : int, optional 559 Sets the number of line segments used to approximate an 560 angle fillet. Note: the use of a `quadsegs` parameter is 561 deprecated and will be gone from the next major release. 562 cap_style : int, optional 563 The styles of caps are: CAP_STYLE.round (1), CAP_STYLE.flat 564 (2), and CAP_STYLE.square (3). 565 join_style : int, optional 566 The styles of joins between offset segments are: 567 JOIN_STYLE.round (1), JOIN_STYLE.mitre (2), and 568 JOIN_STYLE.bevel (3). 569 mitre_limit : float, optional 570 The mitre limit ratio is used for very sharp corners. The 571 mitre ratio is the ratio of the distance from the corner to 572 the end of the mitred offset corner. When two line segments 573 meet at a sharp angle, a miter join will extend the original 574 geometry. To prevent unreasonable geometry, the mitre limit 575 allows controlling the maximum length of the join corner. 576 Corners with a ratio which exceed the limit will be beveled. 577 single_side : bool, optional 578 The side used is determined by the sign of the buffer 579 distance: 580 581 a positive distance indicates the left-hand side 582 a negative distance indicates the right-hand side 583 584 The single-sided buffer of point geometries is the same as 585 the regular buffer. The End Cap Style for single-sided 586 buffers is always ignored, and forced to the equivalent of 587 CAP_FLAT. 588 589 Returns 590 ------- 591 Geometry 592 593 Notes 594 ----- 595 The return value is a strictly two-dimensional geometry. All 596 Z coordinates of the original geometry will be ignored. 597 598 Examples 599 -------- 600 >>> from shapely.wkt import loads 601 >>> g = loads('POINT (0.0 0.0)') 602 >>> g.buffer(1.0).area # 16-gon approx of a unit radius circle 603 3.1365484905459384 604 >>> g.buffer(1.0, 128).area # 128-gon approximation 605 3.141513801144299 606 >>> g.buffer(1.0, 3).area # triangle approximation 607 3.0 608 >>> list(g.buffer(1.0, cap_style=CAP_STYLE.square).exterior.coords) 609 [(1.0, 1.0), (1.0, -1.0), (-1.0, -1.0), (-1.0, 1.0), (1.0, 1.0)] 610 >>> g.buffer(1.0, cap_style=CAP_STYLE.square).area 611 4.0 612 613 """ 614 if quadsegs is not None: 615 warn( 616 "The `quadsegs` argument is deprecated. Use `resolution`.", 617 DeprecationWarning) 618 res = quadsegs 619 else: 620 res = resolution 621 622 if mitre_limit == 0.0: 623 raise ValueError( 624 'Cannot compute offset from zero-length line segment') 625 626 if 'buffer_with_params' in self.impl: 627 params = self._lgeos.GEOSBufferParams_create() 628 self._lgeos.GEOSBufferParams_setEndCapStyle(params, cap_style) 629 self._lgeos.GEOSBufferParams_setJoinStyle(params, join_style) 630 self._lgeos.GEOSBufferParams_setMitreLimit(params, mitre_limit) 631 self._lgeos.GEOSBufferParams_setQuadrantSegments(params, res) 632 self._lgeos.GEOSBufferParams_setSingleSided(params, single_sided) 633 return geom_factory(self.impl['buffer_with_params'](self, params, distance)) 634 635 if cap_style == CAP_STYLE.round and join_style == JOIN_STYLE.round: 636 return geom_factory(self.impl['buffer'](self, distance, res)) 637 638 return geom_factory(self.impl['buffer_with_style'](self, distance, res, 639 cap_style, 640 join_style, 641 mitre_limit)) 642 643 @delegated 644 def simplify(self, tolerance, preserve_topology=True): 645 """Returns a simplified geometry produced by the Douglas-Peucker 646 algorithm 647 648 Coordinates of the simplified geometry will be no more than the 649 tolerance distance from the original. Unless the topology preserving 650 option is used, the algorithm may produce self-intersecting or 651 otherwise invalid geometries. 652 """ 653 if preserve_topology: 654 op = self.impl['topology_preserve_simplify'] 655 else: 656 op = self.impl['simplify'] 657 return geom_factory(op(self, tolerance)) 658 659 def normalize(self): 660 """Converts geometry to normal form (or canonical form). 661 662 This method orders the coordinates, rings of a polygon and parts of 663 multi geometries consistently. Typically useful for testing purposes 664 (for example in combination with `equals_exact`). 665 666 Examples 667 -------- 668 >>> from shapely.wkt import loads 669 >>> p = loads("MULTILINESTRING((0 0, 1 1), (3 3, 2 2))") 670 >>> p.normalize().wkt 671 'MULTILINESTRING ((2 2, 3 3), (0 0, 1 1))' 672 """ 673 # self.impl['normalize'](self) 674 if self._geom is None: 675 raise InvalidGeometryError("Null geometry supports no operations") 676 geom_cloned = lgeos.GEOSGeom_clone(self._geom) 677 lgeos.GEOSNormalize(geom_cloned) 678 return geom_factory(geom_cloned) 679 680 # Binary operations 681 # ----------------- 682 683 def difference(self, other): 684 """Returns the difference of the geometries""" 685 return geom_factory(self.impl['difference'](self, other)) 686 687 def intersection(self, other): 688 """Returns the intersection of the geometries""" 689 return geom_factory(self.impl['intersection'](self, other)) 690 691 def symmetric_difference(self, other): 692 """Returns the symmetric difference of the geometries 693 (Shapely geometry)""" 694 return geom_factory(self.impl['symmetric_difference'](self, other)) 695 696 def union(self, other): 697 """Returns the union of the geometries (Shapely geometry)""" 698 return geom_factory(self.impl['union'](self, other)) 699 700 # Unary predicates 701 # ---------------- 702 703 @property 704 def has_z(self): 705 """True if the geometry's coordinate sequence(s) have z values (are 706 3-dimensional)""" 707 return bool(self.impl['has_z'](self)) 708 709 @property 710 def is_empty(self): 711 """True if the set of points in this geometry is empty, else False""" 712 return (self._geom is None) or bool(self.impl['is_empty'](self)) 713 714 @property 715 def is_ring(self): 716 """True if the geometry is a closed ring, else False""" 717 return bool(self.impl['is_ring'](self)) 718 719 @property 720 def is_closed(self): 721 """True if the geometry is closed, else False 722 723 Applicable only to 1-D geometries.""" 724 if self.geom_type == 'LinearRing': 725 return True 726 elif self.geom_type == 'LineString': 727 if 'is_closed' in self.impl: 728 return bool(self.impl['is_closed'](self)) 729 else: 730 return self.coords[0] == self.coords[-1] 731 else: 732 return False 733 734 @property 735 def is_simple(self): 736 """True if the geometry is simple, meaning that any self-intersections 737 are only at boundary points, else False""" 738 return bool(self.impl['is_simple'](self)) 739 740 @property 741 def is_valid(self): 742 """True if the geometry is valid (definition depends on sub-class), 743 else False""" 744 return bool(self.impl['is_valid'](self)) 745 746 # Binary predicates 747 # ----------------- 748 749 def relate(self, other): 750 """Returns the DE-9IM intersection matrix for the two geometries 751 (string)""" 752 return self.impl['relate'](self, other) 753 754 def covers(self, other): 755 """Returns True if the geometry covers the other, else False""" 756 return bool(self.impl['covers'](self, other)) 757 758 def covered_by(self, other): 759 """Returns True if the geometry is covered by the other, else False""" 760 return bool(self.impl['covered_by'](self, other)) 761 762 def contains(self, other): 763 """Returns True if the geometry contains the other, else False""" 764 return bool(self.impl['contains'](self, other)) 765 766 def crosses(self, other): 767 """Returns True if the geometries cross, else False""" 768 return bool(self.impl['crosses'](self, other)) 769 770 def disjoint(self, other): 771 """Returns True if geometries are disjoint, else False""" 772 return bool(self.impl['disjoint'](self, other)) 773 774 def equals(self, other): 775 """Returns True if geometries are equal, else False. 776 777 This method considers point-set equality (or topological 778 equality), and is equivalent to (self.within(other) & 779 self.contains(other)). 780 781 Examples 782 -------- 783 >>> LineString( 784 ... [(0, 0), (2, 2)] 785 ... ).equals( 786 ... LineString([(0, 0), (1, 1), (2, 2)]) 787 ... ) 788 True 789 790 Returns 791 ------- 792 bool 793 794 """ 795 return bool(self.impl['equals'](self, other)) 796 797 def intersects(self, other): 798 """Returns True if geometries intersect, else False""" 799 return bool(self.impl['intersects'](self, other)) 800 801 def overlaps(self, other): 802 """Returns True if geometries overlap, else False""" 803 return bool(self.impl['overlaps'](self, other)) 804 805 def touches(self, other): 806 """Returns True if geometries touch, else False""" 807 return bool(self.impl['touches'](self, other)) 808 809 def within(self, other): 810 """Returns True if geometry is within the other, else False""" 811 return bool(self.impl['within'](self, other)) 812 813 def equals_exact(self, other, tolerance): 814 """True if geometries are equal to within a specified 815 tolerance. 816 817 Parameters 818 ---------- 819 other : BaseGeometry 820 The other geometry object in this comparison. 821 tolerance : float 822 Absolute tolerance in the same units as coordinates. 823 824 This method considers coordinate equality, which requires 825 coordinates to be equal and in the same order for all components 826 of a geometry. 827 828 Because of this it is possible for "equals()" to be True for two 829 geometries and "equals_exact()" to be False. 830 831 Examples 832 -------- 833 >>> LineString( 834 ... [(0, 0), (2, 2)] 835 ... ).equals_exact( 836 ... LineString([(0, 0), (1, 1), (2, 2)]), 837 ... 1e-6 838 ... ) 839 False 840 841 Returns 842 ------- 843 bool 844 845 """ 846 return bool(self.impl['equals_exact'](self, other, tolerance)) 847 848 def almost_equals(self, other, decimal=6): 849 """True if geometries are equal at all coordinates to a 850 specified decimal place. 851 852 .. deprecated:: 1.8.0 The 'almost_equals()' method is deprecated 853 and will be removed in Shapely 2.0 because the name is 854 confusing. The 'equals_exact()' method should be used 855 instead. 856 857 Refers to approximate coordinate equality, which requires 858 coordinates to be approximately equal and in the same order for 859 all components of a geometry. 860 861 Because of this it is possible for "equals()" to be True for two 862 geometries and "almost_equals()" to be False. 863 864 Examples 865 -------- 866 >>> LineString( 867 ... [(0, 0), (2, 2)] 868 ... ).equals_exact( 869 ... LineString([(0, 0), (1, 1), (2, 2)]), 870 ... 1e-6 871 ... ) 872 False 873 874 Returns 875 ------- 876 bool 877 878 """ 879 warn( 880 "The 'almost_equals()' method is deprecated and will be removed in Shapely 2.0", 881 ShapelyDeprecationWarning, 882 stacklevel=2, 883 ) 884 return self.equals_exact(other, 0.5 * 10 ** (-decimal)) 885 886 def relate_pattern(self, other, pattern): 887 """Returns True if the DE-9IM string code for the relationship between 888 the geometries satisfies the pattern, else False""" 889 pattern = c_char_p(pattern.encode('ascii')) 890 return bool(self.impl['relate_pattern'](self, other, pattern)) 891 892 # Linear referencing 893 # ------------------ 894 895 @delegated 896 def project(self, other, normalized=False): 897 """Returns the distance along this geometry to a point nearest the 898 specified point 899 900 If the normalized arg is True, return the distance normalized to the 901 length of the linear geometry. 902 """ 903 if normalized: 904 op = self.impl['project_normalized'] 905 else: 906 op = self.impl['project'] 907 return op(self, other) 908 909 @delegated 910 @exceptNull 911 def interpolate(self, distance, normalized=False): 912 """Return a point at the specified distance along a linear geometry 913 914 Negative length values are taken as measured in the reverse 915 direction from the end of the geometry. Out-of-range index 916 values are handled by clamping them to the valid range of values. 917 If the normalized arg is True, the distance will be interpreted as a 918 fraction of the geometry's length. 919 """ 920 if normalized: 921 op = self.impl['interpolate_normalized'] 922 else: 923 op = self.impl['interpolate'] 924 return geom_factory(op(self, distance)) 925 926 927class BaseMultipartGeometry(BaseGeometry): 928 929 def shape_factory(self, *args): 930 # Factory for part instances, usually a geometry class 931 raise NotImplementedError("To be implemented by derived classes") 932 933 @property 934 def ctypes(self): 935 raise NotImplementedError( 936 "Multi-part geometries have no ctypes representations") 937 938 @property 939 def __array_interface__(self): 940 """Provide the Numpy array protocol.""" 941 raise AttributeError("Multi-part geometries do not themselves " 942 "provide the array interface") 943 944 def _get_coords(self): 945 raise NotImplementedError("Sub-geometries may have coordinate " 946 "sequences, but collections do not") 947 948 def _set_coords(self, ob): 949 raise NotImplementedError("Sub-geometries may have coordinate " 950 "sequences, but collections do not") 951 952 @property 953 def coords(self): 954 raise NotImplementedError( 955 "Multi-part geometries do not provide a coordinate sequence") 956 957 @property 958 def geoms(self): 959 if self.is_empty: 960 return [] 961 return GeometrySequence(self, self.shape_factory) 962 963 def __bool__(self): 964 return self.is_empty is False 965 966 def __iter__(self): 967 """ 968 .. deprecated:: 1.8 969 Iteration over multi-part geometries is deprecated and will be removed in 970 Shapely 2.0. Use the `geoms` property to access the constituent parts of 971 a multi-part geometry. 972 """ 973 warn( 974 "Iteration over multi-part geometries is deprecated and will be removed in " 975 "Shapely 2.0. Use the `geoms` property to access the constituent parts of " 976 "a multi-part geometry.", ShapelyDeprecationWarning, stacklevel=2) 977 if not self.is_empty: 978 return iter(self.geoms) 979 else: 980 return iter([]) 981 982 def __len__(self): 983 warn( 984 "__len__ for multi-part geometries is deprecated and will be removed in " 985 "Shapely 2.0. Check the length of the `geoms` property instead to get the " 986 " number of parts of a multi-part geometry.", 987 ShapelyDeprecationWarning, stacklevel=2) 988 if not self.is_empty: 989 return len(self.geoms) 990 else: 991 return 0 992 993 def __getitem__(self, index): 994 """ 995 .. deprecated:: 1.8 996 __getitem__ for multi-part geometries is deprecated and will be removed in 997 Shapely 2.0. Use the `geoms` property to access the constituent parts of 998 a multi-part geometry. 999 """ 1000 warn( 1001 "__getitem__ for multi-part geometries is deprecated and will be removed in " 1002 "Shapely 2.0. Use the `geoms` property to access the constituent parts of " 1003 "a multi-part geometry.", ShapelyDeprecationWarning, stacklevel=2) 1004 if not self.is_empty: 1005 return self.geoms[index] 1006 else: 1007 return ()[index] 1008 1009 def __eq__(self, other): 1010 return ( 1011 type(other) == type(self) and 1012 len(self.geoms) == len(other.geoms) and 1013 all(x == y for x, y in zip(self.geoms, other.geoms)) 1014 ) 1015 1016 def __ne__(self, other): 1017 return not self.__eq__(other) 1018 1019 __hash__ = None 1020 1021 def svg(self, scale_factor=1., color=None): 1022 """Returns a group of SVG elements for the multipart geometry. 1023 1024 Parameters 1025 ========== 1026 scale_factor : float 1027 Multiplication factor for the SVG stroke-width. Default is 1. 1028 color : str, optional 1029 Hex string for stroke or fill color. Default is to use "#66cc99" 1030 if geometry is valid, and "#ff3333" if invalid. 1031 """ 1032 if self.is_empty: 1033 return '<g />' 1034 if color is None: 1035 color = "#66cc99" if self.is_valid else "#ff3333" 1036 return '<g>' + \ 1037 ''.join(p.svg(scale_factor, color) for p in self.geoms) + \ 1038 '</g>' 1039 1040 1041class GeometrySequence: 1042 """ 1043 Iterative access to members of a homogeneous multipart geometry. 1044 """ 1045 1046 # Attributes 1047 # ---------- 1048 # _factory : callable 1049 # Returns instances of Shapely geometries 1050 # _geom : c_void_p 1051 # Ctypes pointer to the parent's GEOS geometry 1052 # _ndim : int 1053 # Number of dimensions (2 or 3, generally) 1054 # __p__ : object 1055 # Parent (Shapely) geometry 1056 shape_factory = None 1057 _geom = None 1058 __p__ = None 1059 _ndim = None 1060 1061 def __init__(self, parent, type): 1062 self.shape_factory = type 1063 self.__p__ = parent 1064 1065 def _update(self): 1066 self._geom = self.__p__._geom 1067 self._ndim = self.__p__._ndim 1068 1069 def _get_geom_item(self, i): 1070 g = self.shape_factory() 1071 g._other_owned = True 1072 g._set_geom(lgeos.GEOSGetGeometryN(self._geom, i)) 1073 g._ndim = self._ndim 1074 g.__p__ = self 1075 return g 1076 1077 def __iter__(self): 1078 self._update() 1079 for i in range(self.__len__()): 1080 yield self._get_geom_item(i) 1081 1082 def __len__(self): 1083 self._update() 1084 return lgeos.GEOSGetNumGeometries(self._geom) 1085 1086 def __getitem__(self, key): 1087 self._update() 1088 m = self.__len__() 1089 if isinstance(key, integer_types): 1090 if key + m < 0 or key >= m: 1091 raise IndexError("index out of range") 1092 if key < 0: 1093 i = m + key 1094 else: 1095 i = key 1096 return self._get_geom_item(i) 1097 elif isinstance(key, slice): 1098 if type(self) == HeterogeneousGeometrySequence: 1099 raise GeometryTypeError( 1100 "Heterogenous geometry collections are not sliceable") 1101 res = [] 1102 start, stop, stride = key.indices(m) 1103 for i in range(start, stop, stride): 1104 res.append(self._get_geom_item(i)) 1105 return type(self.__p__)(res or None) 1106 else: 1107 raise TypeError("key must be an index or slice") 1108 1109 @property 1110 def _longest(self): 1111 max = 0 1112 for g in iter(self): 1113 l = len(g.coords) 1114 if l > max: 1115 max = l 1116 1117 1118class HeterogeneousGeometrySequence(GeometrySequence): 1119 """ 1120 Iterative access to a heterogeneous sequence of geometries. 1121 """ 1122 1123 def __init__(self, parent): 1124 super().__init__(parent, None) 1125 1126 def _get_geom_item(self, i): 1127 sub = lgeos.GEOSGetGeometryN(self._geom, i) 1128 g = geom_factory(sub, parent=self) 1129 g._other_owned = True 1130 return g 1131 1132 1133class EmptyGeometry(BaseGeometry): 1134 def __init__(self): 1135 """Create an empty geometry.""" 1136 BaseGeometry.__init__(self) 1137