/django/contrib/gis/geos/geometry.py
Python | 688 lines | 467 code | 68 blank | 153 comment | 47 complexity | 93c15f594a439d9fa5d6c3b7832dd4dd MD5 | raw file
Possible License(s): BSD-3-Clause
1""" 2 This module contains the 'base' GEOSGeometry object -- all GEOS Geometries 3 inherit from this object. 4""" 5# Python, ctypes and types dependencies. 6import re 7import warnings 8from ctypes import addressof, byref, c_double, c_size_t 9 10# super-class for mutable list behavior 11from django.contrib.gis.geos.mutable_list import ListMixin 12 13# GEOS-related dependencies. 14from django.contrib.gis.geos.base import GEOSBase, gdal 15from django.contrib.gis.geos.coordseq import GEOSCoordSeq 16from django.contrib.gis.geos.error import GEOSException, GEOSIndexError 17from django.contrib.gis.geos.libgeos import GEOM_PTR, GEOS_PREPARE 18from django.contrib.gis.geos.mutable_list import ListMixin 19 20# All other functions in this module come from the ctypes 21# prototypes module -- which handles all interaction with 22# the underlying GEOS library. 23from django.contrib.gis.geos import prototypes as capi 24 25# These functions provide access to a thread-local instance 26# of their corresponding GEOS I/O class. 27from django.contrib.gis.geos.prototypes.io import wkt_r, wkt_w, wkb_r, wkb_w, ewkb_w, ewkb_w3d 28 29# For recognizing geometry input. 30from django.contrib.gis.geometry.regex import hex_regex, wkt_regex, json_regex 31 32class GEOSGeometry(GEOSBase, ListMixin): 33 "A class that, generally, encapsulates a GEOS geometry." 34 35 # Raise GEOSIndexError instead of plain IndexError 36 # (see ticket #4740 and GEOSIndexError docstring) 37 _IndexError = GEOSIndexError 38 39 ptr_type = GEOM_PTR 40 41 #### Python 'magic' routines #### 42 def __init__(self, geo_input, srid=None): 43 """ 44 The base constructor for GEOS geometry objects, and may take the 45 following inputs: 46 47 * strings: 48 - WKT 49 - HEXEWKB (a PostGIS-specific canonical form) 50 - GeoJSON (requires GDAL) 51 * buffer: 52 - WKB 53 54 The `srid` keyword is used to specify the Source Reference Identifier 55 (SRID) number for this Geometry. If not set, the SRID will be None. 56 """ 57 if isinstance(geo_input, basestring): 58 if isinstance(geo_input, unicode): 59 # Encoding to ASCII, WKT or HEXEWKB doesn't need any more. 60 geo_input = geo_input.encode('ascii') 61 62 wkt_m = wkt_regex.match(geo_input) 63 if wkt_m: 64 # Handling WKT input. 65 if wkt_m.group('srid'): srid = int(wkt_m.group('srid')) 66 g = wkt_r().read(wkt_m.group('wkt')) 67 elif hex_regex.match(geo_input): 68 # Handling HEXEWKB input. 69 g = wkb_r().read(geo_input) 70 elif gdal.GEOJSON and json_regex.match(geo_input): 71 # Handling GeoJSON input. 72 g = wkb_r().read(gdal.OGRGeometry(geo_input).wkb) 73 else: 74 raise ValueError('String or unicode input unrecognized as WKT EWKT, and HEXEWKB.') 75 elif isinstance(geo_input, GEOM_PTR): 76 # When the input is a pointer to a geomtry (GEOM_PTR). 77 g = geo_input 78 elif isinstance(geo_input, buffer): 79 # When the input is a buffer (WKB). 80 g = wkb_r().read(geo_input) 81 elif isinstance(geo_input, GEOSGeometry): 82 g = capi.geom_clone(geo_input.ptr) 83 else: 84 # Invalid geometry type. 85 raise TypeError('Improper geometry input type: %s' % str(type(geo_input))) 86 87 if bool(g): 88 # Setting the pointer object with a valid pointer. 89 self.ptr = g 90 else: 91 raise GEOSException('Could not initialize GEOS Geometry with given input.') 92 93 # Post-initialization setup. 94 self._post_init(srid) 95 96 def _post_init(self, srid): 97 "Helper routine for performing post-initialization setup." 98 # Setting the SRID, if given. 99 if srid and isinstance(srid, int): self.srid = srid 100 101 # Setting the class type (e.g., Point, Polygon, etc.) 102 self.__class__ = GEOS_CLASSES[self.geom_typeid] 103 104 # Setting the coordinate sequence for the geometry (will be None on 105 # geometries that do not have coordinate sequences) 106 self._set_cs() 107 108 def __del__(self): 109 """ 110 Destroys this Geometry; in other words, frees the memory used by the 111 GEOS C++ object. 112 """ 113 if self._ptr: capi.destroy_geom(self._ptr) 114 115 def __copy__(self): 116 """ 117 Returns a clone because the copy of a GEOSGeometry may contain an 118 invalid pointer location if the original is garbage collected. 119 """ 120 return self.clone() 121 122 def __deepcopy__(self, memodict): 123 """ 124 The `deepcopy` routine is used by the `Node` class of django.utils.tree; 125 thus, the protocol routine needs to be implemented to return correct 126 copies (clones) of these GEOS objects, which use C pointers. 127 """ 128 return self.clone() 129 130 def __str__(self): 131 "WKT is used for the string representation." 132 return self.wkt 133 134 def __repr__(self): 135 "Short-hand representation because WKT may be very large." 136 return '<%s object at %s>' % (self.geom_type, hex(addressof(self.ptr))) 137 138 # Pickling support 139 def __getstate__(self): 140 # The pickled state is simply a tuple of the WKB (in string form) 141 # and the SRID. 142 return str(self.wkb), self.srid 143 144 def __setstate__(self, state): 145 # Instantiating from the tuple state that was pickled. 146 wkb, srid = state 147 ptr = wkb_r().read(buffer(wkb)) 148 if not ptr: raise GEOSException('Invalid Geometry loaded from pickled state.') 149 self.ptr = ptr 150 self._post_init(srid) 151 152 # Comparison operators 153 def __eq__(self, other): 154 """ 155 Equivalence testing, a Geometry may be compared with another Geometry 156 or a WKT representation. 157 """ 158 if isinstance(other, basestring): 159 return self.wkt == other 160 elif isinstance(other, GEOSGeometry): 161 return self.equals_exact(other) 162 else: 163 return False 164 165 def __ne__(self, other): 166 "The not equals operator." 167 return not (self == other) 168 169 ### Geometry set-like operations ### 170 # Thanks to Sean Gillies for inspiration: 171 # http://lists.gispython.org/pipermail/community/2007-July/001034.html 172 # g = g1 | g2 173 def __or__(self, other): 174 "Returns the union of this Geometry and the other." 175 return self.union(other) 176 177 # g = g1 & g2 178 def __and__(self, other): 179 "Returns the intersection of this Geometry and the other." 180 return self.intersection(other) 181 182 # g = g1 - g2 183 def __sub__(self, other): 184 "Return the difference this Geometry and the other." 185 return self.difference(other) 186 187 # g = g1 ^ g2 188 def __xor__(self, other): 189 "Return the symmetric difference of this Geometry and the other." 190 return self.sym_difference(other) 191 192 #### Coordinate Sequence Routines #### 193 @property 194 def has_cs(self): 195 "Returns True if this Geometry has a coordinate sequence, False if not." 196 # Only these geometries are allowed to have coordinate sequences. 197 if isinstance(self, (Point, LineString, LinearRing)): 198 return True 199 else: 200 return False 201 202 def _set_cs(self): 203 "Sets the coordinate sequence for this Geometry." 204 if self.has_cs: 205 self._cs = GEOSCoordSeq(capi.get_cs(self.ptr), self.hasz) 206 else: 207 self._cs = None 208 209 @property 210 def coord_seq(self): 211 "Returns a clone of the coordinate sequence for this Geometry." 212 if self.has_cs: 213 return self._cs.clone() 214 215 #### Geometry Info #### 216 @property 217 def geom_type(self): 218 "Returns a string representing the Geometry type, e.g. 'Polygon'" 219 return capi.geos_type(self.ptr) 220 221 @property 222 def geom_typeid(self): 223 "Returns an integer representing the Geometry type." 224 return capi.geos_typeid(self.ptr) 225 226 @property 227 def num_geom(self): 228 "Returns the number of geometries in the Geometry." 229 return capi.get_num_geoms(self.ptr) 230 231 @property 232 def num_coords(self): 233 "Returns the number of coordinates in the Geometry." 234 return capi.get_num_coords(self.ptr) 235 236 @property 237 def num_points(self): 238 "Returns the number points, or coordinates, in the Geometry." 239 return self.num_coords 240 241 @property 242 def dims(self): 243 "Returns the dimension of this Geometry (0=point, 1=line, 2=surface)." 244 return capi.get_dims(self.ptr) 245 246 def normalize(self): 247 "Converts this Geometry to normal form (or canonical form)." 248 return capi.geos_normalize(self.ptr) 249 250 #### Unary predicates #### 251 @property 252 def empty(self): 253 """ 254 Returns a boolean indicating whether the set of points in this Geometry 255 are empty. 256 """ 257 return capi.geos_isempty(self.ptr) 258 259 @property 260 def hasz(self): 261 "Returns whether the geometry has a 3D dimension." 262 return capi.geos_hasz(self.ptr) 263 264 @property 265 def ring(self): 266 "Returns whether or not the geometry is a ring." 267 return capi.geos_isring(self.ptr) 268 269 @property 270 def simple(self): 271 "Returns false if the Geometry not simple." 272 return capi.geos_issimple(self.ptr) 273 274 @property 275 def valid(self): 276 "This property tests the validity of this Geometry." 277 return capi.geos_isvalid(self.ptr) 278 279 @property 280 def valid_reason(self): 281 """ 282 Returns a string containing the reason for any invalidity. 283 """ 284 if not GEOS_PREPARE: 285 raise GEOSException('Upgrade GEOS to 3.1 to get validity reason.') 286 return capi.geos_isvalidreason(self.ptr) 287 288 #### Binary predicates. #### 289 def contains(self, other): 290 "Returns true if other.within(this) returns true." 291 return capi.geos_contains(self.ptr, other.ptr) 292 293 def crosses(self, other): 294 """ 295 Returns true if the DE-9IM intersection matrix for the two Geometries 296 is T*T****** (for a point and a curve,a point and an area or a line and 297 an area) 0******** (for two curves). 298 """ 299 return capi.geos_crosses(self.ptr, other.ptr) 300 301 def disjoint(self, other): 302 """ 303 Returns true if the DE-9IM intersection matrix for the two Geometries 304 is FF*FF****. 305 """ 306 return capi.geos_disjoint(self.ptr, other.ptr) 307 308 def equals(self, other): 309 """ 310 Returns true if the DE-9IM intersection matrix for the two Geometries 311 is T*F**FFF*. 312 """ 313 return capi.geos_equals(self.ptr, other.ptr) 314 315 def equals_exact(self, other, tolerance=0): 316 """ 317 Returns true if the two Geometries are exactly equal, up to a 318 specified tolerance. 319 """ 320 return capi.geos_equalsexact(self.ptr, other.ptr, float(tolerance)) 321 322 def intersects(self, other): 323 "Returns true if disjoint returns false." 324 return capi.geos_intersects(self.ptr, other.ptr) 325 326 def overlaps(self, other): 327 """ 328 Returns true if the DE-9IM intersection matrix for the two Geometries 329 is T*T***T** (for two points or two surfaces) 1*T***T** (for two curves). 330 """ 331 return capi.geos_overlaps(self.ptr, other.ptr) 332 333 def relate_pattern(self, other, pattern): 334 """ 335 Returns true if the elements in the DE-9IM intersection matrix for the 336 two Geometries match the elements in pattern. 337 """ 338 if not isinstance(pattern, basestring) or len(pattern) > 9: 339 raise GEOSException('invalid intersection matrix pattern') 340 return capi.geos_relatepattern(self.ptr, other.ptr, pattern) 341 342 def touches(self, other): 343 """ 344 Returns true if the DE-9IM intersection matrix for the two Geometries 345 is FT*******, F**T***** or F***T****. 346 """ 347 return capi.geos_touches(self.ptr, other.ptr) 348 349 def within(self, other): 350 """ 351 Returns true if the DE-9IM intersection matrix for the two Geometries 352 is T*F**F***. 353 """ 354 return capi.geos_within(self.ptr, other.ptr) 355 356 #### SRID Routines #### 357 def get_srid(self): 358 "Gets the SRID for the geometry, returns None if no SRID is set." 359 s = capi.geos_get_srid(self.ptr) 360 if s == 0: return None 361 else: return s 362 363 def set_srid(self, srid): 364 "Sets the SRID for the geometry." 365 capi.geos_set_srid(self.ptr, srid) 366 srid = property(get_srid, set_srid) 367 368 #### Output Routines #### 369 @property 370 def ewkt(self): 371 """ 372 Returns the EWKT (WKT + SRID) of the Geometry. Note that Z values 373 are *not* included in this representation because GEOS does not yet 374 support serializing them. 375 """ 376 if self.get_srid(): return 'SRID=%s;%s' % (self.srid, self.wkt) 377 else: return self.wkt 378 379 @property 380 def wkt(self): 381 "Returns the WKT (Well-Known Text) representation of this Geometry." 382 return wkt_w().write(self) 383 384 @property 385 def hex(self): 386 """ 387 Returns the WKB of this Geometry in hexadecimal form. Please note 388 that the SRID and Z values are not included in this representation 389 because it is not a part of the OGC specification (use the `hexewkb` 390 property instead). 391 """ 392 # A possible faster, all-python, implementation: 393 # str(self.wkb).encode('hex') 394 return wkb_w().write_hex(self) 395 396 @property 397 def hexewkb(self): 398 """ 399 Returns the EWKB of this Geometry in hexadecimal form. This is an 400 extension of the WKB specification that includes SRID and Z values 401 that are a part of this geometry. 402 """ 403 if self.hasz: 404 if not GEOS_PREPARE: 405 # See: http://trac.osgeo.org/geos/ticket/216 406 raise GEOSException('Upgrade GEOS to 3.1 to get valid 3D HEXEWKB.') 407 return ewkb_w3d().write_hex(self) 408 else: 409 return ewkb_w().write_hex(self) 410 411 @property 412 def json(self): 413 """ 414 Returns GeoJSON representation of this Geometry if GDAL 1.5+ 415 is installed. 416 """ 417 if gdal.GEOJSON: 418 return self.ogr.json 419 else: 420 raise GEOSException('GeoJSON output only supported on GDAL 1.5+.') 421 geojson = json 422 423 @property 424 def wkb(self): 425 """ 426 Returns the WKB (Well-Known Binary) representation of this Geometry 427 as a Python buffer. SRID and Z values are not included, use the 428 `ewkb` property instead. 429 """ 430 return wkb_w().write(self) 431 432 @property 433 def ewkb(self): 434 """ 435 Return the EWKB representation of this Geometry as a Python buffer. 436 This is an extension of the WKB specification that includes any SRID 437 and Z values that are a part of this geometry. 438 """ 439 if self.hasz: 440 if not GEOS_PREPARE: 441 # See: http://trac.osgeo.org/geos/ticket/216 442 raise GEOSException('Upgrade GEOS to 3.1 to get valid 3D EWKB.') 443 return ewkb_w3d().write(self) 444 else: 445 return ewkb_w().write(self) 446 447 @property 448 def kml(self): 449 "Returns the KML representation of this Geometry." 450 gtype = self.geom_type 451 return '<%s>%s</%s>' % (gtype, self.coord_seq.kml, gtype) 452 453 @property 454 def prepared(self): 455 """ 456 Returns a PreparedGeometry corresponding to this geometry -- it is 457 optimized for the contains, intersects, and covers operations. 458 """ 459 if GEOS_PREPARE: 460 return PreparedGeometry(self) 461 else: 462 raise GEOSException('GEOS 3.1+ required for prepared geometry support.') 463 464 #### GDAL-specific output routines #### 465 @property 466 def ogr(self): 467 "Returns the OGR Geometry for this Geometry." 468 if gdal.HAS_GDAL: 469 if self.srid: 470 return gdal.OGRGeometry(self.wkb, self.srid) 471 else: 472 return gdal.OGRGeometry(self.wkb) 473 else: 474 raise GEOSException('GDAL required to convert to an OGRGeometry.') 475 476 @property 477 def srs(self): 478 "Returns the OSR SpatialReference for SRID of this Geometry." 479 if gdal.HAS_GDAL: 480 if self.srid: 481 return gdal.SpatialReference(self.srid) 482 else: 483 return None 484 else: 485 raise GEOSException('GDAL required to return a SpatialReference object.') 486 487 @property 488 def crs(self): 489 "Alias for `srs` property." 490 return self.srs 491 492 def transform(self, ct, clone=False): 493 """ 494 Requires GDAL. Transforms the geometry according to the given 495 transformation object, which may be an integer SRID, and WKT or 496 PROJ.4 string. By default, the geometry is transformed in-place and 497 nothing is returned. However if the `clone` keyword is set, then this 498 geometry will not be modified and a transformed clone will be returned 499 instead. 500 """ 501 srid = self.srid 502 503 if ct == srid: 504 # short-circuit where source & dest SRIDs match 505 if clone: 506 return self.clone() 507 else: 508 return 509 510 if (srid is None) or (srid < 0): 511 warnings.warn("Calling transform() with no SRID set does no transformation!", 512 stacklevel=2) 513 warnings.warn("Calling transform() with no SRID will raise GEOSException in v1.5", 514 FutureWarning, stacklevel=2) 515 return 516 517 if not gdal.HAS_GDAL: 518 raise GEOSException("GDAL library is not available to transform() geometry.") 519 520 # Creating an OGR Geometry, which is then transformed. 521 g = gdal.OGRGeometry(self.wkb, srid) 522 g.transform(ct) 523 # Getting a new GEOS pointer 524 ptr = wkb_r().read(g.wkb) 525 if clone: 526 # User wants a cloned transformed geometry returned. 527 return GEOSGeometry(ptr, srid=g.srid) 528 if ptr: 529 # Reassigning pointer, and performing post-initialization setup 530 # again due to the reassignment. 531 capi.destroy_geom(self.ptr) 532 self.ptr = ptr 533 self._post_init(g.srid) 534 else: 535 raise GEOSException('Transformed WKB was invalid.') 536 537 #### Topology Routines #### 538 def _topology(self, gptr): 539 "Helper routine to return Geometry from the given pointer." 540 return GEOSGeometry(gptr, srid=self.srid) 541 542 @property 543 def boundary(self): 544 "Returns the boundary as a newly allocated Geometry object." 545 return self._topology(capi.geos_boundary(self.ptr)) 546 547 def buffer(self, width, quadsegs=8): 548 """ 549 Returns a geometry that represents all points whose distance from this 550 Geometry is less than or equal to distance. Calculations are in the 551 Spatial Reference System of this Geometry. The optional third parameter sets 552 the number of segment used to approximate a quarter circle (defaults to 8). 553 (Text from PostGIS documentation at ch. 6.1.3) 554 """ 555 return self._topology(capi.geos_buffer(self.ptr, width, quadsegs)) 556 557 @property 558 def centroid(self): 559 """ 560 The centroid is equal to the centroid of the set of component Geometries 561 of highest dimension (since the lower-dimension geometries contribute zero 562 "weight" to the centroid). 563 """ 564 return self._topology(capi.geos_centroid(self.ptr)) 565 566 @property 567 def convex_hull(self): 568 """ 569 Returns the smallest convex Polygon that contains all the points 570 in the Geometry. 571 """ 572 return self._topology(capi.geos_convexhull(self.ptr)) 573 574 def difference(self, other): 575 """ 576 Returns a Geometry representing the points making up this Geometry 577 that do not make up other. 578 """ 579 return self._topology(capi.geos_difference(self.ptr, other.ptr)) 580 581 @property 582 def envelope(self): 583 "Return the envelope for this geometry (a polygon)." 584 return self._topology(capi.geos_envelope(self.ptr)) 585 586 def intersection(self, other): 587 "Returns a Geometry representing the points shared by this Geometry and other." 588 return self._topology(capi.geos_intersection(self.ptr, other.ptr)) 589 590 @property 591 def point_on_surface(self): 592 "Computes an interior point of this Geometry." 593 return self._topology(capi.geos_pointonsurface(self.ptr)) 594 595 def relate(self, other): 596 "Returns the DE-9IM intersection matrix for this Geometry and the other." 597 return capi.geos_relate(self.ptr, other.ptr) 598 599 def simplify(self, tolerance=0.0, preserve_topology=False): 600 """ 601 Returns the Geometry, simplified using the Douglas-Peucker algorithm 602 to the specified tolerance (higher tolerance => less points). If no 603 tolerance provided, defaults to 0. 604 605 By default, this function does not preserve topology - e.g. polygons can 606 be split, collapse to lines or disappear holes can be created or 607 disappear, and lines can cross. By specifying preserve_topology=True, 608 the result will have the same dimension and number of components as the 609 input. This is significantly slower. 610 """ 611 if preserve_topology: 612 return self._topology(capi.geos_preservesimplify(self.ptr, tolerance)) 613 else: 614 return self._topology(capi.geos_simplify(self.ptr, tolerance)) 615 616 def sym_difference(self, other): 617 """ 618 Returns a set combining the points in this Geometry not in other, 619 and the points in other not in this Geometry. 620 """ 621 return self._topology(capi.geos_symdifference(self.ptr, other.ptr)) 622 623 def union(self, other): 624 "Returns a Geometry representing all the points in this Geometry and other." 625 return self._topology(capi.geos_union(self.ptr, other.ptr)) 626 627 #### Other Routines #### 628 @property 629 def area(self): 630 "Returns the area of the Geometry." 631 return capi.geos_area(self.ptr, byref(c_double())) 632 633 def distance(self, other): 634 """ 635 Returns the distance between the closest points on this Geometry 636 and the other. Units will be in those of the coordinate system of 637 the Geometry. 638 """ 639 if not isinstance(other, GEOSGeometry): 640 raise TypeError('distance() works only on other GEOS Geometries.') 641 return capi.geos_distance(self.ptr, other.ptr, byref(c_double())) 642 643 @property 644 def extent(self): 645 """ 646 Returns the extent of this geometry as a 4-tuple, consisting of 647 (xmin, ymin, xmax, ymax). 648 """ 649 env = self.envelope 650 if isinstance(env, Point): 651 xmin, ymin = env.tuple 652 xmax, ymax = xmin, ymin 653 else: 654 xmin, ymin = env[0][0] 655 xmax, ymax = env[0][2] 656 return (xmin, ymin, xmax, ymax) 657 658 @property 659 def length(self): 660 """ 661 Returns the length of this Geometry (e.g., 0 for point, or the 662 circumfrence of a Polygon). 663 """ 664 return capi.geos_length(self.ptr, byref(c_double())) 665 666 def clone(self): 667 "Clones this Geometry." 668 return GEOSGeometry(capi.geom_clone(self.ptr), srid=self.srid) 669 670# Class mapping dictionary. Has to be at the end to avoid import 671# conflicts with GEOSGeometry. 672from django.contrib.gis.geos.linestring import LineString, LinearRing 673from django.contrib.gis.geos.point import Point 674from django.contrib.gis.geos.polygon import Polygon 675from django.contrib.gis.geos.collections import GeometryCollection, MultiPoint, MultiLineString, MultiPolygon 676GEOS_CLASSES = {0 : Point, 677 1 : LineString, 678 2 : LinearRing, 679 3 : Polygon, 680 4 : MultiPoint, 681 5 : MultiLineString, 682 6 : MultiPolygon, 683 7 : GeometryCollection, 684 } 685 686# If supported, import the PreparedGeometry class. 687if GEOS_PREPARE: 688 from django.contrib.gis.geos.prepared import PreparedGeometry