You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
643 lines
22 KiB
643 lines
22 KiB
""" |
|
The OGRGeometry is a wrapper for using the OGR Geometry class |
|
(see http://www.gdal.org/ogr/classOGRGeometry.html). OGRGeometry |
|
may be instantiated when reading geometries from OGR Data Sources |
|
(e.g. SHP files), or when given OGC WKT (a string). |
|
|
|
While the 'full' API is not present yet, the API is "pythonic" unlike |
|
the traditional and "next-generation" OGR Python bindings. One major |
|
advantage OGR Geometries have over their GEOS counterparts is support |
|
for spatial reference systems and their transformation. |
|
|
|
Example: |
|
>>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference |
|
>>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)' |
|
>>> pnt = OGRGeometry(wkt1) |
|
>>> print pnt |
|
POINT (-90 30) |
|
>>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84')) |
|
>>> mpnt.add(wkt1) |
|
>>> mpnt.add(wkt1) |
|
>>> print mpnt |
|
MULTIPOINT (-90 30,-90 30) |
|
>>> print mpnt.srs.name |
|
WGS 84 |
|
>>> print mpnt.srs.proj |
|
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs |
|
>>> mpnt.transform_to(SpatialReference('NAD27')) |
|
>>> print mpnt.proj |
|
+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs |
|
>>> print mpnt |
|
MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641) |
|
|
|
The OGRGeomType class is to make it easy to specify an OGR geometry type: |
|
>>> from django.contrib.gis.gdal import OGRGeomType |
|
>>> gt1 = OGRGeomType(3) # Using an integer for the type |
|
>>> gt2 = OGRGeomType('Polygon') # Using a string |
|
>>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive |
|
>>> print gt1 == 3, gt1 == 'Polygon' # Equivalence works w/non-OGRGeomType objects |
|
True |
|
""" |
|
# Python library requisites. |
|
import re, sys |
|
from binascii import a2b_hex |
|
from ctypes import byref, string_at, c_char_p, c_double, c_ubyte, c_void_p |
|
from types import UnicodeType |
|
|
|
# Getting GDAL prerequisites |
|
from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope |
|
from django.contrib.gis.gdal.error import OGRException, OGRIndexError, SRSException |
|
from django.contrib.gis.gdal.geomtype import OGRGeomType |
|
from django.contrib.gis.gdal.srs import SpatialReference, CoordTransform |
|
|
|
# Getting the ctypes prototype functions that interface w/the GDAL C library. |
|
from django.contrib.gis.gdal.prototypes.geom import * |
|
from django.contrib.gis.gdal.prototypes.srs import clone_srs |
|
|
|
# For more information, see the OGR C API source code: |
|
# http://www.gdal.org/ogr/ogr__api_8h.html |
|
# |
|
# The OGR_G_* routines are relevant here. |
|
|
|
# Regular expressions for recognizing HEXEWKB and WKT. |
|
hex_regex = re.compile(r'^[0-9A-F]+$', re.I) |
|
wkt_regex = re.compile(r'^(?P<type>POINT|LINESTRING|LINEARRING|POLYGON|MULTIPOINT|MULTILINESTRING|MULTIPOLYGON|GEOMETRYCOLLECTION)[ACEGIMLONPSRUTY\d,\.\-\(\) ]+$', re.I) |
|
json_regex = re.compile(r'^\{[\s\w,\-\.\"\'\:\[\]]+\}$') |
|
|
|
#### OGRGeometry Class #### |
|
class OGRGeometry(object): |
|
"Generally encapsulates an OGR geometry." |
|
|
|
def __init__(self, geom_input, srs=None): |
|
"Initializes Geometry on either WKT or an OGR pointer as input." |
|
|
|
self._ptr = c_void_p(None) # Initially NULL |
|
str_instance = isinstance(geom_input, basestring) |
|
|
|
# If HEX, unpack input to to a binary buffer. |
|
if str_instance and hex_regex.match(geom_input): |
|
geom_input = buffer(a2b_hex(geom_input.upper())) |
|
str_instance = False |
|
|
|
# Constructing the geometry, |
|
if str_instance: |
|
# Checking if unicode |
|
if isinstance(geom_input, UnicodeType): |
|
# Encoding to ASCII, WKT or HEX doesn't need any more. |
|
geo_input = geo_input.encode('ascii') |
|
|
|
wkt_m = wkt_regex.match(geom_input) |
|
json_m = json_regex.match(geom_input) |
|
if wkt_m: |
|
if wkt_m.group('type').upper() == 'LINEARRING': |
|
# OGR_G_CreateFromWkt doesn't work with LINEARRING WKT. |
|
# See http://trac.osgeo.org/gdal/ticket/1992. |
|
g = create_geom(OGRGeomType(wkt_m.group('type')).num) |
|
import_wkt(g, byref(c_char_p(geom_input))) |
|
else: |
|
g = from_wkt(byref(c_char_p(geom_input)), None, byref(c_void_p())) |
|
elif json_m: |
|
if GEOJSON: |
|
g = from_json(geom_input) |
|
else: |
|
raise NotImplementedError('GeoJSON input only supported on GDAL 1.5+.') |
|
else: |
|
# Seeing if the input is a valid short-hand string |
|
# (e.g., 'Point', 'POLYGON'). |
|
ogr_t = OGRGeomType(geom_input) |
|
g = create_geom(OGRGeomType(geom_input).num) |
|
elif isinstance(geom_input, buffer): |
|
# WKB was passed in |
|
g = from_wkb(str(geom_input), None, byref(c_void_p()), len(geom_input)) |
|
elif isinstance(geom_input, OGRGeomType): |
|
# OGRGeomType was passed in, an empty geometry will be created. |
|
g = create_geom(geom_input.num) |
|
elif isinstance(geom_input, c_void_p): |
|
# OGR pointer (c_void_p) was the input. |
|
g = geom_input |
|
else: |
|
raise OGRException('Invalid input type for OGR Geometry construction: %s' % type(geom_input)) |
|
|
|
# Now checking the Geometry pointer before finishing initialization |
|
# by setting the pointer for the object. |
|
if not g: |
|
raise OGRException('Cannot create OGR Geometry from input: %s' % str(geom_input)) |
|
self._ptr = g |
|
|
|
# Assigning the SpatialReference object to the geometry, if valid. |
|
if bool(srs): self.srs = srs |
|
|
|
# Setting the class depending upon the OGR Geometry Type |
|
self.__class__ = GEO_CLASSES[self.geom_type.num] |
|
|
|
def __del__(self): |
|
"Deletes this Geometry." |
|
if self._ptr: destroy_geom(self._ptr) |
|
|
|
### Geometry set-like operations ### |
|
# g = g1 | g2 |
|
def __or__(self, other): |
|
"Returns the union of the two geometries." |
|
return self.union(other) |
|
|
|
# g = g1 & g2 |
|
def __and__(self, other): |
|
"Returns the intersection of this Geometry and the other." |
|
return self.intersection(other) |
|
|
|
# g = g1 - g2 |
|
def __sub__(self, other): |
|
"Return the difference this Geometry and the other." |
|
return self.difference(other) |
|
|
|
# g = g1 ^ g2 |
|
def __xor__(self, other): |
|
"Return the symmetric difference of this Geometry and the other." |
|
return self.sym_difference(other) |
|
|
|
def __eq__(self, other): |
|
"Is this Geometry equal to the other?" |
|
return self.equals(other) |
|
|
|
def __ne__(self, other): |
|
"Tests for inequality." |
|
return not self.equals(other) |
|
|
|
def __str__(self): |
|
"WKT is used for the string representation." |
|
return self.wkt |
|
|
|
#### Geometry Properties #### |
|
@property |
|
def dimension(self): |
|
"Returns 0 for points, 1 for lines, and 2 for surfaces." |
|
return get_dims(self._ptr) |
|
|
|
@property |
|
def coord_dim(self): |
|
"Returns the coordinate dimension of the Geometry." |
|
return get_coord_dims(self._ptr) |
|
|
|
@property |
|
def geom_count(self): |
|
"The number of elements in this Geometry." |
|
return get_geom_count(self._ptr) |
|
|
|
@property |
|
def point_count(self): |
|
"Returns the number of Points in this Geometry." |
|
return get_point_count(self._ptr) |
|
|
|
@property |
|
def num_points(self): |
|
"Alias for `point_count` (same name method in GEOS API.)" |
|
return self.point_count |
|
|
|
@property |
|
def num_coords(self): |
|
"Alais for `point_count`." |
|
return self.point_count |
|
|
|
@property |
|
def geom_type(self): |
|
"Returns the Type for this Geometry." |
|
try: |
|
return OGRGeomType(get_geom_type(self._ptr)) |
|
except OGRException: |
|
# VRT datasources return an invalid geometry type |
|
# number, but a valid name -- we'll try that instead. |
|
# See: http://trac.osgeo.org/gdal/ticket/2491 |
|
return OGRGeomType(get_geom_name(self._ptr)) |
|
|
|
@property |
|
def geom_name(self): |
|
"Returns the Name of this Geometry." |
|
return get_geom_name(self._ptr) |
|
|
|
@property |
|
def area(self): |
|
"Returns the area for a LinearRing, Polygon, or MultiPolygon; 0 otherwise." |
|
return get_area(self._ptr) |
|
|
|
@property |
|
def envelope(self): |
|
"Returns the envelope for this Geometry." |
|
# TODO: Fix Envelope() for Point geometries. |
|
return Envelope(get_envelope(self._ptr, byref(OGREnvelope()))) |
|
|
|
@property |
|
def extent(self): |
|
"Returns the envelope as a 4-tuple, instead of as an Envelope object." |
|
return self.envelope.tuple |
|
|
|
#### SpatialReference-related Properties #### |
|
|
|
# The SRS property |
|
def get_srs(self): |
|
"Returns the Spatial Reference for this Geometry." |
|
try: |
|
srs_ptr = get_geom_srs(self._ptr) |
|
return SpatialReference(clone_srs(srs_ptr)) |
|
except SRSException: |
|
return None |
|
|
|
def set_srs(self, srs): |
|
"Sets the SpatialReference for this geometry." |
|
if isinstance(srs, SpatialReference): |
|
srs_ptr = clone_srs(srs._ptr) |
|
elif isinstance(srs, (int, long, basestring)): |
|
sr = SpatialReference(srs) |
|
srs_ptr = clone_srs(sr._ptr) |
|
else: |
|
raise TypeError('Cannot assign spatial reference with object of type: %s' % type(srs)) |
|
assign_srs(self._ptr, srs_ptr) |
|
|
|
srs = property(get_srs, set_srs) |
|
|
|
# The SRID property |
|
def get_srid(self): |
|
if self.srs: return self.srs.srid |
|
else: return None |
|
|
|
def set_srid(self, srid): |
|
if isinstance(srid, (int, long)): |
|
self.srs = srid |
|
else: |
|
raise TypeError('SRID must be set with an integer.') |
|
|
|
srid = property(get_srid, set_srid) |
|
|
|
#### Output Methods #### |
|
@property |
|
def geos(self): |
|
"Returns a GEOSGeometry object from this OGRGeometry." |
|
from django.contrib.gis.geos import GEOSGeometry |
|
return GEOSGeometry(self.wkb, self.srid) |
|
|
|
@property |
|
def gml(self): |
|
"Returns the GML representation of the Geometry." |
|
return to_gml(self._ptr) |
|
|
|
@property |
|
def hex(self): |
|
"Returns the hexadecimal representation of the WKB (a string)." |
|
return str(self.wkb).encode('hex').upper() |
|
#return b2a_hex(self.wkb).upper() |
|
|
|
@property |
|
def json(self): |
|
if GEOJSON: |
|
return to_json(self._ptr) |
|
else: |
|
raise NotImplementedError('GeoJSON output only supported on GDAL 1.5+.') |
|
geojson = json |
|
|
|
@property |
|
def wkb_size(self): |
|
"Returns the size of the WKB buffer." |
|
return get_wkbsize(self._ptr) |
|
|
|
@property |
|
def wkb(self): |
|
"Returns the WKB representation of the Geometry." |
|
if sys.byteorder == 'little': |
|
byteorder = 1 # wkbNDR (from ogr_core.h) |
|
else: |
|
byteorder = 0 # wkbXDR |
|
sz = self.wkb_size |
|
# Creating the unsigned character buffer, and passing it in by reference. |
|
buf = (c_ubyte * sz)() |
|
wkb = to_wkb(self._ptr, byteorder, byref(buf)) |
|
# Returning a buffer of the string at the pointer. |
|
return buffer(string_at(buf, sz)) |
|
|
|
@property |
|
def wkt(self): |
|
"Returns the WKT representation of the Geometry." |
|
return to_wkt(self._ptr, byref(c_char_p())) |
|
|
|
#### Geometry Methods #### |
|
def clone(self): |
|
"Clones this OGR Geometry." |
|
return OGRGeometry(clone_geom(self._ptr), self.srs) |
|
|
|
def close_rings(self): |
|
""" |
|
If there are any rings within this geometry that have not been |
|
closed, this routine will do so by adding the starting point at the |
|
end. |
|
""" |
|
# Closing the open rings. |
|
geom_close_rings(self._ptr) |
|
|
|
def transform(self, coord_trans, clone=False): |
|
""" |
|
Transforms this geometry to a different spatial reference system. |
|
May take a CoordTransform object, a SpatialReference object, string |
|
WKT or PROJ.4, and/or an integer SRID. By default nothing is returned |
|
and the geometry is transformed in-place. However, if the `clone` |
|
keyword is set, then a transformed clone of this geometry will be |
|
returned. |
|
""" |
|
if clone: |
|
klone = self.clone() |
|
klone.transform(coord_trans) |
|
return klone |
|
if isinstance(coord_trans, CoordTransform): |
|
geom_transform(self._ptr, coord_trans._ptr) |
|
elif isinstance(coord_trans, SpatialReference): |
|
geom_transform_to(self._ptr, coord_trans._ptr) |
|
elif isinstance(coord_trans, (int, long, basestring)): |
|
sr = SpatialReference(coord_trans) |
|
geom_transform_to(self._ptr, sr._ptr) |
|
else: |
|
raise TypeError('Transform only accepts CoordTransform, SpatialReference, string, and integer objects.') |
|
|
|
def transform_to(self, srs): |
|
"For backwards-compatibility." |
|
self.transform(srs) |
|
|
|
#### Topology Methods #### |
|
def _topology(self, func, other): |
|
"""A generalized function for topology operations, takes a GDAL function and |
|
the other geometry to perform the operation on.""" |
|
if not isinstance(other, OGRGeometry): |
|
raise TypeError('Must use another OGRGeometry object for topology operations!') |
|
|
|
# Returning the output of the given function with the other geometry's |
|
# pointer. |
|
return func(self._ptr, other._ptr) |
|
|
|
def intersects(self, other): |
|
"Returns True if this geometry intersects with the other." |
|
return self._topology(ogr_intersects, other) |
|
|
|
def equals(self, other): |
|
"Returns True if this geometry is equivalent to the other." |
|
return self._topology(ogr_equals, other) |
|
|
|
def disjoint(self, other): |
|
"Returns True if this geometry and the other are spatially disjoint." |
|
return self._topology(ogr_disjoint, other) |
|
|
|
def touches(self, other): |
|
"Returns True if this geometry touches the other." |
|
return self._topology(ogr_touches, other) |
|
|
|
def crosses(self, other): |
|
"Returns True if this geometry crosses the other." |
|
return self._topology(ogr_crosses, other) |
|
|
|
def within(self, other): |
|
"Returns True if this geometry is within the other." |
|
return self._topology(ogr_within, other) |
|
|
|
def contains(self, other): |
|
"Returns True if this geometry contains the other." |
|
return self._topology(ogr_contains, other) |
|
|
|
def overlaps(self, other): |
|
"Returns True if this geometry overlaps the other." |
|
return self._topology(ogr_overlaps, other) |
|
|
|
#### Geometry-generation Methods #### |
|
def _geomgen(self, gen_func, other=None): |
|
"A helper routine for the OGR routines that generate geometries." |
|
if isinstance(other, OGRGeometry): |
|
return OGRGeometry(gen_func(self._ptr, other._ptr), self.srs) |
|
else: |
|
return OGRGeometry(gen_func(self._ptr), self.srs) |
|
|
|
@property |
|
def boundary(self): |
|
"Returns the boundary of this geometry." |
|
return self._geomgen(get_boundary) |
|
|
|
@property |
|
def convex_hull(self): |
|
""" |
|
Returns the smallest convex Polygon that contains all the points in |
|
this Geometry. |
|
""" |
|
return self._geomgen(geom_convex_hull) |
|
|
|
def difference(self, other): |
|
""" |
|
Returns a new geometry consisting of the region which is the difference |
|
of this geometry and the other. |
|
""" |
|
return self._geomgen(geom_diff, other) |
|
|
|
def intersection(self, other): |
|
""" |
|
Returns a new geometry consisting of the region of intersection of this |
|
geometry and the other. |
|
""" |
|
return self._geomgen(geom_intersection, other) |
|
|
|
def sym_difference(self, other): |
|
""" |
|
Returns a new geometry which is the symmetric difference of this |
|
geometry and the other. |
|
""" |
|
return self._geomgen(geom_sym_diff, other) |
|
|
|
def union(self, other): |
|
""" |
|
Returns a new geometry consisting of the region which is the union of |
|
this geometry and the other. |
|
""" |
|
return self._geomgen(geom_union, other) |
|
|
|
# The subclasses for OGR Geometry. |
|
class Point(OGRGeometry): |
|
|
|
@property |
|
def x(self): |
|
"Returns the X coordinate for this Point." |
|
return getx(self._ptr, 0) |
|
|
|
@property |
|
def y(self): |
|
"Returns the Y coordinate for this Point." |
|
return gety(self._ptr, 0) |
|
|
|
@property |
|
def z(self): |
|
"Returns the Z coordinate for this Point." |
|
if self.coord_dim == 3: |
|
return getz(self._ptr, 0) |
|
|
|
@property |
|
def tuple(self): |
|
"Returns the tuple of this point." |
|
if self.coord_dim == 2: |
|
return (self.x, self.y) |
|
elif self.coord_dim == 3: |
|
return (self.x, self.y, self.z) |
|
coords = tuple |
|
|
|
class LineString(OGRGeometry): |
|
|
|
def __getitem__(self, index): |
|
"Returns the Point at the given index." |
|
if index >= 0 and index < self.point_count: |
|
x, y, z = c_double(), c_double(), c_double() |
|
get_point(self._ptr, index, byref(x), byref(y), byref(z)) |
|
dim = self.coord_dim |
|
if dim == 1: |
|
return (x.value,) |
|
elif dim == 2: |
|
return (x.value, y.value) |
|
elif dim == 3: |
|
return (x.value, y.value, z.value) |
|
else: |
|
raise OGRIndexError('index out of range: %s' % str(index)) |
|
|
|
def __iter__(self): |
|
"Iterates over each point in the LineString." |
|
for i in xrange(self.point_count): |
|
yield self[i] |
|
|
|
def __len__(self): |
|
"The length returns the number of points in the LineString." |
|
return self.point_count |
|
|
|
@property |
|
def tuple(self): |
|
"Returns the tuple representation of this LineString." |
|
return tuple([self[i] for i in xrange(len(self))]) |
|
coords = tuple |
|
|
|
def _listarr(self, func): |
|
""" |
|
Internal routine that returns a sequence (list) corresponding with |
|
the given function. |
|
""" |
|
return [func(self._ptr, i) for i in xrange(len(self))] |
|
|
|
@property |
|
def x(self): |
|
"Returns the X coordinates in a list." |
|
return self._listarr(getx) |
|
|
|
@property |
|
def y(self): |
|
"Returns the Y coordinates in a list." |
|
return self._listarr(gety) |
|
|
|
@property |
|
def z(self): |
|
"Returns the Z coordinates in a list." |
|
if self.coord_dim == 3: |
|
return self._listarr(getz) |
|
|
|
# LinearRings are used in Polygons. |
|
class LinearRing(LineString): pass |
|
|
|
class Polygon(OGRGeometry): |
|
|
|
def __len__(self): |
|
"The number of interior rings in this Polygon." |
|
return self.geom_count |
|
|
|
def __iter__(self): |
|
"Iterates through each ring in the Polygon." |
|
for i in xrange(self.geom_count): |
|
yield self[i] |
|
|
|
def __getitem__(self, index): |
|
"Gets the ring at the specified index." |
|
if index < 0 or index >= self.geom_count: |
|
raise OGRIndexError('index out of range: %s' % index) |
|
else: |
|
return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs) |
|
|
|
# Polygon Properties |
|
@property |
|
def shell(self): |
|
"Returns the shell of this Polygon." |
|
return self[0] # First ring is the shell |
|
exterior_ring = shell |
|
|
|
@property |
|
def tuple(self): |
|
"Returns a tuple of LinearRing coordinate tuples." |
|
return tuple([self[i].tuple for i in xrange(self.geom_count)]) |
|
coords = tuple |
|
|
|
@property |
|
def point_count(self): |
|
"The number of Points in this Polygon." |
|
# Summing up the number of points in each ring of the Polygon. |
|
return sum([self[i].point_count for i in xrange(self.geom_count)]) |
|
|
|
@property |
|
def centroid(self): |
|
"Returns the centroid (a Point) of this Polygon." |
|
# The centroid is a Point, create a geometry for this. |
|
p = OGRGeometry(OGRGeomType('Point')) |
|
get_centroid(self._ptr, p._ptr) |
|
return p |
|
|
|
# Geometry Collection base class. |
|
class GeometryCollection(OGRGeometry): |
|
"The Geometry Collection class." |
|
|
|
def __getitem__(self, index): |
|
"Gets the Geometry at the specified index." |
|
if index < 0 or index >= self.geom_count: |
|
raise OGRIndexError('index out of range: %s' % index) |
|
else: |
|
return OGRGeometry(clone_geom(get_geom_ref(self._ptr, index)), self.srs) |
|
|
|
def __iter__(self): |
|
"Iterates over each Geometry." |
|
for i in xrange(self.geom_count): |
|
yield self[i] |
|
|
|
def __len__(self): |
|
"The number of geometries in this Geometry Collection." |
|
return self.geom_count |
|
|
|
def add(self, geom): |
|
"Add the geometry to this Geometry Collection." |
|
if isinstance(geom, OGRGeometry): |
|
if isinstance(geom, self.__class__): |
|
for g in geom: add_geom(self._ptr, g._ptr) |
|
else: |
|
add_geom(self._ptr, geom._ptr) |
|
elif isinstance(geom, basestring): |
|
tmp = OGRGeometry(geom) |
|
add_geom(self._ptr, tmp._ptr) |
|
else: |
|
raise OGRException('Must add an OGRGeometry.') |
|
|
|
@property |
|
def point_count(self): |
|
"The number of Points in this Geometry Collection." |
|
# Summing up the number of points in each geometry in this collection |
|
return sum([self[i].point_count for i in xrange(self.geom_count)]) |
|
|
|
@property |
|
def tuple(self): |
|
"Returns a tuple representation of this Geometry Collection." |
|
return tuple([self[i].tuple for i in xrange(self.geom_count)]) |
|
coords = tuple |
|
|
|
# Multiple Geometry types. |
|
class MultiPoint(GeometryCollection): pass |
|
class MultiLineString(GeometryCollection): pass |
|
class MultiPolygon(GeometryCollection): pass |
|
|
|
# Class mapping dictionary (using the OGRwkbGeometryType as the key) |
|
GEO_CLASSES = {1 : Point, |
|
2 : LineString, |
|
3 : Polygon, |
|
4 : MultiPoint, |
|
5 : MultiLineString, |
|
6 : MultiPolygon, |
|
7 : GeometryCollection, |
|
101: LinearRing, |
|
}
|
|
|