|
0
|
1 |
""" |
|
|
2 |
The OGRGeometry is a wrapper for using the OGR Geometry class |
|
|
3 |
(see http://www.gdal.org/ogr/classOGRGeometry.html). OGRGeometry |
|
|
4 |
may be instantiated when reading geometries from OGR Data Sources |
|
|
5 |
(e.g. SHP files), or when given OGC WKT (a string). |
|
|
6 |
|
|
|
7 |
While the 'full' API is not present yet, the API is "pythonic" unlike |
|
|
8 |
the traditional and "next-generation" OGR Python bindings. One major |
|
|
9 |
advantage OGR Geometries have over their GEOS counterparts is support |
|
|
10 |
for spatial reference systems and their transformation. |
|
|
11 |
|
|
|
12 |
Example: |
|
|
13 |
>>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference |
|
|
14 |
>>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)' |
|
|
15 |
>>> pnt = OGRGeometry(wkt1) |
|
|
16 |
>>> print pnt |
|
|
17 |
POINT (-90 30) |
|
|
18 |
>>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84')) |
|
|
19 |
>>> mpnt.add(wkt1) |
|
|
20 |
>>> mpnt.add(wkt1) |
|
|
21 |
>>> print mpnt |
|
|
22 |
MULTIPOINT (-90 30,-90 30) |
|
|
23 |
>>> print mpnt.srs.name |
|
|
24 |
WGS 84 |
|
|
25 |
>>> print mpnt.srs.proj |
|
|
26 |
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs |
|
|
27 |
>>> mpnt.transform_to(SpatialReference('NAD27')) |
|
|
28 |
>>> print mpnt.proj |
|
|
29 |
+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs |
|
|
30 |
>>> print mpnt |
|
|
31 |
MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641) |
|
29
|
32 |
|
|
0
|
33 |
The OGRGeomType class is to make it easy to specify an OGR geometry type: |
|
|
34 |
>>> from django.contrib.gis.gdal import OGRGeomType |
|
|
35 |
>>> gt1 = OGRGeomType(3) # Using an integer for the type |
|
|
36 |
>>> gt2 = OGRGeomType('Polygon') # Using a string |
|
|
37 |
>>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive |
|
|
38 |
>>> print gt1 == 3, gt1 == 'Polygon' # Equivalence works w/non-OGRGeomType objects |
|
|
39 |
True |
|
|
40 |
""" |
|
|
41 |
# Python library requisites. |
|
29
|
42 |
import sys |
|
0
|
43 |
from binascii import a2b_hex |
|
|
44 |
from ctypes import byref, string_at, c_char_p, c_double, c_ubyte, c_void_p |
|
|
45 |
|
|
|
46 |
# Getting GDAL prerequisites |
|
|
47 |
from django.contrib.gis.gdal.base import GDALBase |
|
|
48 |
from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope |
|
|
49 |
from django.contrib.gis.gdal.error import OGRException, OGRIndexError, SRSException |
|
|
50 |
from django.contrib.gis.gdal.geomtype import OGRGeomType |
|
29
|
51 |
from django.contrib.gis.gdal.libgdal import GEOJSON, GDAL_VERSION |
|
0
|
52 |
from django.contrib.gis.gdal.srs import SpatialReference, CoordTransform |
|
|
53 |
|
|
|
54 |
# Getting the ctypes prototype functions that interface w/the GDAL C library. |
|
|
55 |
from django.contrib.gis.gdal.prototypes import geom as capi, srs as srs_api |
|
29
|
56 |
|
|
|
57 |
# For recognizing geometry input. |
|
|
58 |
from django.contrib.gis.geometry.regex import hex_regex, wkt_regex, json_regex |
|
0
|
59 |
|
|
|
60 |
# For more information, see the OGR C API source code: |
|
|
61 |
# http://www.gdal.org/ogr/ogr__api_8h.html |
|
|
62 |
# |
|
|
63 |
# The OGR_G_* routines are relevant here. |
|
|
64 |
|
|
|
65 |
#### OGRGeometry Class #### |
|
|
66 |
class OGRGeometry(GDALBase): |
|
|
67 |
"Generally encapsulates an OGR geometry." |
|
|
68 |
|
|
|
69 |
def __init__(self, geom_input, srs=None): |
|
|
70 |
"Initializes Geometry on either WKT or an OGR pointer as input." |
|
|
71 |
|
|
|
72 |
str_instance = isinstance(geom_input, basestring) |
|
|
73 |
|
|
|
74 |
# If HEX, unpack input to to a binary buffer. |
|
|
75 |
if str_instance and hex_regex.match(geom_input): |
|
|
76 |
geom_input = buffer(a2b_hex(geom_input.upper())) |
|
|
77 |
str_instance = False |
|
|
78 |
|
|
29
|
79 |
# Constructing the geometry, |
|
0
|
80 |
if str_instance: |
|
|
81 |
# Checking if unicode |
|
|
82 |
if isinstance(geom_input, unicode): |
|
|
83 |
# Encoding to ASCII, WKT or HEX doesn't need any more. |
|
|
84 |
geom_input = geom_input.encode('ascii') |
|
|
85 |
|
|
|
86 |
wkt_m = wkt_regex.match(geom_input) |
|
|
87 |
json_m = json_regex.match(geom_input) |
|
|
88 |
if wkt_m: |
|
29
|
89 |
if wkt_m.group('srid'): |
|
|
90 |
# If there's EWKT, set the SRS w/value of the SRID. |
|
|
91 |
srs = int(wkt_m.group('srid')) |
|
0
|
92 |
if wkt_m.group('type').upper() == 'LINEARRING': |
|
|
93 |
# OGR_G_CreateFromWkt doesn't work with LINEARRING WKT. |
|
|
94 |
# See http://trac.osgeo.org/gdal/ticket/1992. |
|
|
95 |
g = capi.create_geom(OGRGeomType(wkt_m.group('type')).num) |
|
29
|
96 |
capi.import_wkt(g, byref(c_char_p(wkt_m.group('wkt')))) |
|
0
|
97 |
else: |
|
29
|
98 |
g = capi.from_wkt(byref(c_char_p(wkt_m.group('wkt'))), None, byref(c_void_p())) |
|
0
|
99 |
elif json_m: |
|
|
100 |
if GEOJSON: |
|
|
101 |
g = capi.from_json(geom_input) |
|
|
102 |
else: |
|
|
103 |
raise NotImplementedError('GeoJSON input only supported on GDAL 1.5+.') |
|
|
104 |
else: |
|
|
105 |
# Seeing if the input is a valid short-hand string |
|
|
106 |
# (e.g., 'Point', 'POLYGON'). |
|
|
107 |
ogr_t = OGRGeomType(geom_input) |
|
|
108 |
g = capi.create_geom(OGRGeomType(geom_input).num) |
|
|
109 |
elif isinstance(geom_input, buffer): |
|
|
110 |
# WKB was passed in |
|
|
111 |
g = capi.from_wkb(str(geom_input), None, byref(c_void_p()), len(geom_input)) |
|
|
112 |
elif isinstance(geom_input, OGRGeomType): |
|
|
113 |
# OGRGeomType was passed in, an empty geometry will be created. |
|
|
114 |
g = capi.create_geom(geom_input.num) |
|
|
115 |
elif isinstance(geom_input, self.ptr_type): |
|
|
116 |
# OGR pointer (c_void_p) was the input. |
|
|
117 |
g = geom_input |
|
|
118 |
else: |
|
|
119 |
raise OGRException('Invalid input type for OGR Geometry construction: %s' % type(geom_input)) |
|
|
120 |
|
|
|
121 |
# Now checking the Geometry pointer before finishing initialization |
|
|
122 |
# by setting the pointer for the object. |
|
|
123 |
if not g: |
|
|
124 |
raise OGRException('Cannot create OGR Geometry from input: %s' % str(geom_input)) |
|
|
125 |
self.ptr = g |
|
|
126 |
|
|
|
127 |
# Assigning the SpatialReference object to the geometry, if valid. |
|
|
128 |
if bool(srs): self.srs = srs |
|
|
129 |
|
|
|
130 |
# Setting the class depending upon the OGR Geometry Type |
|
|
131 |
self.__class__ = GEO_CLASSES[self.geom_type.num] |
|
|
132 |
|
|
29
|
133 |
def __del__(self): |
|
|
134 |
"Deletes this Geometry." |
|
|
135 |
if self._ptr: capi.destroy_geom(self._ptr) |
|
|
136 |
|
|
|
137 |
# Pickle routines |
|
|
138 |
def __getstate__(self): |
|
|
139 |
srs = self.srs |
|
|
140 |
if srs: |
|
|
141 |
srs = srs.wkt |
|
|
142 |
else: |
|
|
143 |
srs = None |
|
|
144 |
return str(self.wkb), srs |
|
|
145 |
|
|
|
146 |
def __setstate__(self, state): |
|
|
147 |
wkb, srs = state |
|
|
148 |
ptr = capi.from_wkb(wkb, None, byref(c_void_p()), len(wkb)) |
|
|
149 |
if not ptr: raise OGRException('Invalid OGRGeometry loaded from pickled state.') |
|
|
150 |
self.ptr = ptr |
|
|
151 |
self.srs = srs |
|
|
152 |
|
|
0
|
153 |
@classmethod |
|
29
|
154 |
def from_bbox(cls, bbox): |
|
0
|
155 |
"Constructs a Polygon from a bounding box (4-tuple)." |
|
|
156 |
x0, y0, x1, y1 = bbox |
|
|
157 |
return OGRGeometry( 'POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))' % ( |
|
|
158 |
x0, y0, x0, y1, x1, y1, x1, y0, x0, y0) ) |
|
|
159 |
|
|
|
160 |
### Geometry set-like operations ### |
|
|
161 |
# g = g1 | g2 |
|
|
162 |
def __or__(self, other): |
|
|
163 |
"Returns the union of the two geometries." |
|
|
164 |
return self.union(other) |
|
|
165 |
|
|
|
166 |
# g = g1 & g2 |
|
|
167 |
def __and__(self, other): |
|
|
168 |
"Returns the intersection of this Geometry and the other." |
|
|
169 |
return self.intersection(other) |
|
|
170 |
|
|
|
171 |
# g = g1 - g2 |
|
|
172 |
def __sub__(self, other): |
|
|
173 |
"Return the difference this Geometry and the other." |
|
|
174 |
return self.difference(other) |
|
|
175 |
|
|
|
176 |
# g = g1 ^ g2 |
|
|
177 |
def __xor__(self, other): |
|
|
178 |
"Return the symmetric difference of this Geometry and the other." |
|
|
179 |
return self.sym_difference(other) |
|
|
180 |
|
|
|
181 |
def __eq__(self, other): |
|
|
182 |
"Is this Geometry equal to the other?" |
|
29
|
183 |
if isinstance(other, OGRGeometry): |
|
|
184 |
return self.equals(other) |
|
|
185 |
else: |
|
|
186 |
return False |
|
0
|
187 |
|
|
|
188 |
def __ne__(self, other): |
|
|
189 |
"Tests for inequality." |
|
29
|
190 |
return not (self == other) |
|
0
|
191 |
|
|
|
192 |
def __str__(self): |
|
|
193 |
"WKT is used for the string representation." |
|
|
194 |
return self.wkt |
|
|
195 |
|
|
|
196 |
#### Geometry Properties #### |
|
|
197 |
@property |
|
|
198 |
def dimension(self): |
|
|
199 |
"Returns 0 for points, 1 for lines, and 2 for surfaces." |
|
|
200 |
return capi.get_dims(self.ptr) |
|
|
201 |
|
|
29
|
202 |
def _get_coord_dim(self): |
|
0
|
203 |
"Returns the coordinate dimension of the Geometry." |
|
29
|
204 |
if isinstance(self, GeometryCollection) and GDAL_VERSION < (1, 5, 2): |
|
|
205 |
# On GDAL versions prior to 1.5.2, there exists a bug in which |
|
|
206 |
# the coordinate dimension of geometry collections is always 2: |
|
|
207 |
# http://trac.osgeo.org/gdal/ticket/2334 |
|
|
208 |
# Here we workaround by returning the coordinate dimension of the |
|
|
209 |
# first geometry in the collection instead. |
|
|
210 |
if len(self): |
|
|
211 |
return capi.get_coord_dim(capi.get_geom_ref(self.ptr, 0)) |
|
|
212 |
return capi.get_coord_dim(self.ptr) |
|
|
213 |
|
|
|
214 |
def _set_coord_dim(self, dim): |
|
|
215 |
"Sets the coordinate dimension of this Geometry." |
|
|
216 |
if not dim in (2, 3): |
|
|
217 |
raise ValueError('Geometry dimension must be either 2 or 3') |
|
|
218 |
capi.set_coord_dim(self.ptr, dim) |
|
|
219 |
|
|
|
220 |
coord_dim = property(_get_coord_dim, _set_coord_dim) |
|
0
|
221 |
|
|
|
222 |
@property |
|
|
223 |
def geom_count(self): |
|
|
224 |
"The number of elements in this Geometry." |
|
|
225 |
return capi.get_geom_count(self.ptr) |
|
|
226 |
|
|
|
227 |
@property |
|
|
228 |
def point_count(self): |
|
|
229 |
"Returns the number of Points in this Geometry." |
|
|
230 |
return capi.get_point_count(self.ptr) |
|
|
231 |
|
|
|
232 |
@property |
|
|
233 |
def num_points(self): |
|
|
234 |
"Alias for `point_count` (same name method in GEOS API.)" |
|
|
235 |
return self.point_count |
|
|
236 |
|
|
|
237 |
@property |
|
|
238 |
def num_coords(self): |
|
|
239 |
"Alais for `point_count`." |
|
|
240 |
return self.point_count |
|
|
241 |
|
|
|
242 |
@property |
|
|
243 |
def geom_type(self): |
|
|
244 |
"Returns the Type for this Geometry." |
|
29
|
245 |
return OGRGeomType(capi.get_geom_type(self.ptr)) |
|
0
|
246 |
|
|
|
247 |
@property |
|
|
248 |
def geom_name(self): |
|
|
249 |
"Returns the Name of this Geometry." |
|
|
250 |
return capi.get_geom_name(self.ptr) |
|
|
251 |
|
|
|
252 |
@property |
|
|
253 |
def area(self): |
|
|
254 |
"Returns the area for a LinearRing, Polygon, or MultiPolygon; 0 otherwise." |
|
|
255 |
return capi.get_area(self.ptr) |
|
|
256 |
|
|
|
257 |
@property |
|
|
258 |
def envelope(self): |
|
|
259 |
"Returns the envelope for this Geometry." |
|
|
260 |
# TODO: Fix Envelope() for Point geometries. |
|
|
261 |
return Envelope(capi.get_envelope(self.ptr, byref(OGREnvelope()))) |
|
|
262 |
|
|
|
263 |
@property |
|
|
264 |
def extent(self): |
|
|
265 |
"Returns the envelope as a 4-tuple, instead of as an Envelope object." |
|
|
266 |
return self.envelope.tuple |
|
|
267 |
|
|
|
268 |
#### SpatialReference-related Properties #### |
|
29
|
269 |
|
|
0
|
270 |
# The SRS property |
|
|
271 |
def _get_srs(self): |
|
|
272 |
"Returns the Spatial Reference for this Geometry." |
|
|
273 |
try: |
|
|
274 |
srs_ptr = capi.get_geom_srs(self.ptr) |
|
|
275 |
return SpatialReference(srs_api.clone_srs(srs_ptr)) |
|
|
276 |
except SRSException: |
|
|
277 |
return None |
|
|
278 |
|
|
|
279 |
def _set_srs(self, srs): |
|
|
280 |
"Sets the SpatialReference for this geometry." |
|
29
|
281 |
# Do not have to clone the `SpatialReference` object pointer because |
|
|
282 |
# when it is assigned to this `OGRGeometry` it's internal OGR |
|
|
283 |
# reference count is incremented, and will likewise be released |
|
|
284 |
# (decremented) when this geometry's destructor is called. |
|
0
|
285 |
if isinstance(srs, SpatialReference): |
|
29
|
286 |
srs_ptr = srs.ptr |
|
0
|
287 |
elif isinstance(srs, (int, long, basestring)): |
|
|
288 |
sr = SpatialReference(srs) |
|
29
|
289 |
srs_ptr = sr.ptr |
|
0
|
290 |
else: |
|
|
291 |
raise TypeError('Cannot assign spatial reference with object of type: %s' % type(srs)) |
|
|
292 |
capi.assign_srs(self.ptr, srs_ptr) |
|
|
293 |
|
|
|
294 |
srs = property(_get_srs, _set_srs) |
|
|
295 |
|
|
|
296 |
# The SRID property |
|
|
297 |
def _get_srid(self): |
|
|
298 |
srs = self.srs |
|
|
299 |
if srs: return srs.srid |
|
|
300 |
return None |
|
|
301 |
|
|
|
302 |
def _set_srid(self, srid): |
|
|
303 |
if isinstance(srid, (int, long)): |
|
|
304 |
self.srs = srid |
|
|
305 |
else: |
|
|
306 |
raise TypeError('SRID must be set with an integer.') |
|
|
307 |
|
|
|
308 |
srid = property(_get_srid, _set_srid) |
|
|
309 |
|
|
|
310 |
#### Output Methods #### |
|
|
311 |
@property |
|
|
312 |
def geos(self): |
|
|
313 |
"Returns a GEOSGeometry object from this OGRGeometry." |
|
|
314 |
from django.contrib.gis.geos import GEOSGeometry |
|
|
315 |
return GEOSGeometry(self.wkb, self.srid) |
|
|
316 |
|
|
|
317 |
@property |
|
|
318 |
def gml(self): |
|
|
319 |
"Returns the GML representation of the Geometry." |
|
|
320 |
return capi.to_gml(self.ptr) |
|
|
321 |
|
|
|
322 |
@property |
|
|
323 |
def hex(self): |
|
|
324 |
"Returns the hexadecimal representation of the WKB (a string)." |
|
|
325 |
return str(self.wkb).encode('hex').upper() |
|
|
326 |
#return b2a_hex(self.wkb).upper() |
|
|
327 |
|
|
|
328 |
@property |
|
|
329 |
def json(self): |
|
|
330 |
""" |
|
|
331 |
Returns the GeoJSON representation of this Geometry (requires |
|
|
332 |
GDAL 1.5+). |
|
|
333 |
""" |
|
29
|
334 |
if GEOJSON: |
|
0
|
335 |
return capi.to_json(self.ptr) |
|
|
336 |
else: |
|
|
337 |
raise NotImplementedError('GeoJSON output only supported on GDAL 1.5+.') |
|
|
338 |
geojson = json |
|
|
339 |
|
|
|
340 |
@property |
|
|
341 |
def kml(self): |
|
|
342 |
"Returns the KML representation of the Geometry." |
|
|
343 |
if GEOJSON: |
|
|
344 |
return capi.to_kml(self.ptr, None) |
|
|
345 |
else: |
|
|
346 |
raise NotImplementedError('KML output only supported on GDAL 1.5+.') |
|
|
347 |
|
|
|
348 |
@property |
|
|
349 |
def wkb_size(self): |
|
|
350 |
"Returns the size of the WKB buffer." |
|
|
351 |
return capi.get_wkbsize(self.ptr) |
|
|
352 |
|
|
|
353 |
@property |
|
|
354 |
def wkb(self): |
|
|
355 |
"Returns the WKB representation of the Geometry." |
|
|
356 |
if sys.byteorder == 'little': |
|
|
357 |
byteorder = 1 # wkbNDR (from ogr_core.h) |
|
|
358 |
else: |
|
|
359 |
byteorder = 0 # wkbXDR |
|
|
360 |
sz = self.wkb_size |
|
|
361 |
# Creating the unsigned character buffer, and passing it in by reference. |
|
|
362 |
buf = (c_ubyte * sz)() |
|
|
363 |
wkb = capi.to_wkb(self.ptr, byteorder, byref(buf)) |
|
|
364 |
# Returning a buffer of the string at the pointer. |
|
|
365 |
return buffer(string_at(buf, sz)) |
|
|
366 |
|
|
|
367 |
@property |
|
|
368 |
def wkt(self): |
|
|
369 |
"Returns the WKT representation of the Geometry." |
|
|
370 |
return capi.to_wkt(self.ptr, byref(c_char_p())) |
|
29
|
371 |
|
|
|
372 |
@property |
|
|
373 |
def ewkt(self): |
|
|
374 |
"Returns the EWKT representation of the Geometry." |
|
|
375 |
srs = self.srs |
|
|
376 |
if srs and srs.srid: |
|
|
377 |
return 'SRID=%s;%s' % (srs.srid, self.wkt) |
|
|
378 |
else: |
|
|
379 |
return self.wkt |
|
|
380 |
|
|
0
|
381 |
#### Geometry Methods #### |
|
|
382 |
def clone(self): |
|
|
383 |
"Clones this OGR Geometry." |
|
|
384 |
return OGRGeometry(capi.clone_geom(self.ptr), self.srs) |
|
|
385 |
|
|
|
386 |
def close_rings(self): |
|
|
387 |
""" |
|
|
388 |
If there are any rings within this geometry that have not been |
|
|
389 |
closed, this routine will do so by adding the starting point at the |
|
|
390 |
end. |
|
|
391 |
""" |
|
|
392 |
# Closing the open rings. |
|
|
393 |
capi.geom_close_rings(self.ptr) |
|
|
394 |
|
|
|
395 |
def transform(self, coord_trans, clone=False): |
|
|
396 |
""" |
|
|
397 |
Transforms this geometry to a different spatial reference system. |
|
|
398 |
May take a CoordTransform object, a SpatialReference object, string |
|
|
399 |
WKT or PROJ.4, and/or an integer SRID. By default nothing is returned |
|
|
400 |
and the geometry is transformed in-place. However, if the `clone` |
|
|
401 |
keyword is set, then a transformed clone of this geometry will be |
|
|
402 |
returned. |
|
|
403 |
""" |
|
|
404 |
if clone: |
|
|
405 |
klone = self.clone() |
|
|
406 |
klone.transform(coord_trans) |
|
|
407 |
return klone |
|
29
|
408 |
|
|
|
409 |
# Have to get the coordinate dimension of the original geometry |
|
|
410 |
# so it can be used to reset the transformed geometry's dimension |
|
|
411 |
# afterwards. This is done because of GDAL bug (in versions prior |
|
|
412 |
# to 1.7) that turns geometries 3D after transformation, see: |
|
|
413 |
# http://trac.osgeo.org/gdal/changeset/17792 |
|
|
414 |
if GDAL_VERSION < (1, 7): |
|
|
415 |
orig_dim = self.coord_dim |
|
|
416 |
|
|
|
417 |
# Depending on the input type, use the appropriate OGR routine |
|
|
418 |
# to perform the transformation. |
|
0
|
419 |
if isinstance(coord_trans, CoordTransform): |
|
|
420 |
capi.geom_transform(self.ptr, coord_trans.ptr) |
|
|
421 |
elif isinstance(coord_trans, SpatialReference): |
|
|
422 |
capi.geom_transform_to(self.ptr, coord_trans.ptr) |
|
|
423 |
elif isinstance(coord_trans, (int, long, basestring)): |
|
|
424 |
sr = SpatialReference(coord_trans) |
|
|
425 |
capi.geom_transform_to(self.ptr, sr.ptr) |
|
|
426 |
else: |
|
29
|
427 |
raise TypeError('Transform only accepts CoordTransform, ' |
|
|
428 |
'SpatialReference, string, and integer objects.') |
|
|
429 |
|
|
|
430 |
# Setting with original dimension, see comment above. |
|
|
431 |
if GDAL_VERSION < (1, 7): |
|
|
432 |
if isinstance(self, GeometryCollection): |
|
|
433 |
# With geometry collections have to set dimension on |
|
|
434 |
# each internal geometry reference, as the collection |
|
|
435 |
# dimension isn't affected. |
|
|
436 |
for i in xrange(len(self)): |
|
|
437 |
internal_ptr = capi.get_geom_ref(self.ptr, i) |
|
|
438 |
if orig_dim != capi.get_coord_dim(internal_ptr): |
|
|
439 |
capi.set_coord_dim(internal_ptr, orig_dim) |
|
|
440 |
else: |
|
|
441 |
if self.coord_dim != orig_dim: |
|
|
442 |
self.coord_dim = orig_dim |
|
0
|
443 |
|
|
|
444 |
def transform_to(self, srs): |
|
|
445 |
"For backwards-compatibility." |
|
|
446 |
self.transform(srs) |
|
|
447 |
|
|
|
448 |
#### Topology Methods #### |
|
|
449 |
def _topology(self, func, other): |
|
|
450 |
"""A generalized function for topology operations, takes a GDAL function and |
|
|
451 |
the other geometry to perform the operation on.""" |
|
|
452 |
if not isinstance(other, OGRGeometry): |
|
|
453 |
raise TypeError('Must use another OGRGeometry object for topology operations!') |
|
|
454 |
|
|
|
455 |
# Returning the output of the given function with the other geometry's |
|
|
456 |
# pointer. |
|
|
457 |
return func(self.ptr, other.ptr) |
|
|
458 |
|
|
|
459 |
def intersects(self, other): |
|
|
460 |
"Returns True if this geometry intersects with the other." |
|
|
461 |
return self._topology(capi.ogr_intersects, other) |
|
29
|
462 |
|
|
0
|
463 |
def equals(self, other): |
|
|
464 |
"Returns True if this geometry is equivalent to the other." |
|
|
465 |
return self._topology(capi.ogr_equals, other) |
|
|
466 |
|
|
|
467 |
def disjoint(self, other): |
|
|
468 |
"Returns True if this geometry and the other are spatially disjoint." |
|
|
469 |
return self._topology(capi.ogr_disjoint, other) |
|
|
470 |
|
|
|
471 |
def touches(self, other): |
|
|
472 |
"Returns True if this geometry touches the other." |
|
|
473 |
return self._topology(capi.ogr_touches, other) |
|
|
474 |
|
|
|
475 |
def crosses(self, other): |
|
|
476 |
"Returns True if this geometry crosses the other." |
|
|
477 |
return self._topology(capi.ogr_crosses, other) |
|
|
478 |
|
|
|
479 |
def within(self, other): |
|
|
480 |
"Returns True if this geometry is within the other." |
|
|
481 |
return self._topology(capi.ogr_within, other) |
|
|
482 |
|
|
|
483 |
def contains(self, other): |
|
|
484 |
"Returns True if this geometry contains the other." |
|
|
485 |
return self._topology(capi.ogr_contains, other) |
|
|
486 |
|
|
|
487 |
def overlaps(self, other): |
|
|
488 |
"Returns True if this geometry overlaps the other." |
|
|
489 |
return self._topology(capi.ogr_overlaps, other) |
|
|
490 |
|
|
|
491 |
#### Geometry-generation Methods #### |
|
|
492 |
def _geomgen(self, gen_func, other=None): |
|
|
493 |
"A helper routine for the OGR routines that generate geometries." |
|
|
494 |
if isinstance(other, OGRGeometry): |
|
|
495 |
return OGRGeometry(gen_func(self.ptr, other.ptr), self.srs) |
|
|
496 |
else: |
|
|
497 |
return OGRGeometry(gen_func(self.ptr), self.srs) |
|
|
498 |
|
|
|
499 |
@property |
|
|
500 |
def boundary(self): |
|
|
501 |
"Returns the boundary of this geometry." |
|
|
502 |
return self._geomgen(capi.get_boundary) |
|
|
503 |
|
|
|
504 |
@property |
|
|
505 |
def convex_hull(self): |
|
|
506 |
""" |
|
29
|
507 |
Returns the smallest convex Polygon that contains all the points in |
|
0
|
508 |
this Geometry. |
|
|
509 |
""" |
|
|
510 |
return self._geomgen(capi.geom_convex_hull) |
|
|
511 |
|
|
|
512 |
def difference(self, other): |
|
|
513 |
""" |
|
|
514 |
Returns a new geometry consisting of the region which is the difference |
|
|
515 |
of this geometry and the other. |
|
|
516 |
""" |
|
|
517 |
return self._geomgen(capi.geom_diff, other) |
|
|
518 |
|
|
|
519 |
def intersection(self, other): |
|
|
520 |
""" |
|
|
521 |
Returns a new geometry consisting of the region of intersection of this |
|
|
522 |
geometry and the other. |
|
|
523 |
""" |
|
|
524 |
return self._geomgen(capi.geom_intersection, other) |
|
|
525 |
|
|
|
526 |
def sym_difference(self, other): |
|
29
|
527 |
""" |
|
0
|
528 |
Returns a new geometry which is the symmetric difference of this |
|
|
529 |
geometry and the other. |
|
|
530 |
""" |
|
|
531 |
return self._geomgen(capi.geom_sym_diff, other) |
|
|
532 |
|
|
|
533 |
def union(self, other): |
|
|
534 |
""" |
|
|
535 |
Returns a new geometry consisting of the region which is the union of |
|
|
536 |
this geometry and the other. |
|
|
537 |
""" |
|
|
538 |
return self._geomgen(capi.geom_union, other) |
|
|
539 |
|
|
|
540 |
# The subclasses for OGR Geometry. |
|
|
541 |
class Point(OGRGeometry): |
|
|
542 |
|
|
|
543 |
@property |
|
|
544 |
def x(self): |
|
|
545 |
"Returns the X coordinate for this Point." |
|
|
546 |
return capi.getx(self.ptr, 0) |
|
|
547 |
|
|
|
548 |
@property |
|
|
549 |
def y(self): |
|
|
550 |
"Returns the Y coordinate for this Point." |
|
|
551 |
return capi.gety(self.ptr, 0) |
|
|
552 |
|
|
|
553 |
@property |
|
|
554 |
def z(self): |
|
|
555 |
"Returns the Z coordinate for this Point." |
|
|
556 |
if self.coord_dim == 3: |
|
|
557 |
return capi.getz(self.ptr, 0) |
|
|
558 |
|
|
|
559 |
@property |
|
|
560 |
def tuple(self): |
|
|
561 |
"Returns the tuple of this point." |
|
|
562 |
if self.coord_dim == 2: |
|
|
563 |
return (self.x, self.y) |
|
|
564 |
elif self.coord_dim == 3: |
|
|
565 |
return (self.x, self.y, self.z) |
|
|
566 |
coords = tuple |
|
|
567 |
|
|
|
568 |
class LineString(OGRGeometry): |
|
|
569 |
|
|
|
570 |
def __getitem__(self, index): |
|
|
571 |
"Returns the Point at the given index." |
|
|
572 |
if index >= 0 and index < self.point_count: |
|
|
573 |
x, y, z = c_double(), c_double(), c_double() |
|
|
574 |
capi.get_point(self.ptr, index, byref(x), byref(y), byref(z)) |
|
|
575 |
dim = self.coord_dim |
|
|
576 |
if dim == 1: |
|
|
577 |
return (x.value,) |
|
|
578 |
elif dim == 2: |
|
|
579 |
return (x.value, y.value) |
|
|
580 |
elif dim == 3: |
|
|
581 |
return (x.value, y.value, z.value) |
|
|
582 |
else: |
|
|
583 |
raise OGRIndexError('index out of range: %s' % str(index)) |
|
|
584 |
|
|
|
585 |
def __iter__(self): |
|
|
586 |
"Iterates over each point in the LineString." |
|
|
587 |
for i in xrange(self.point_count): |
|
|
588 |
yield self[i] |
|
|
589 |
|
|
|
590 |
def __len__(self): |
|
|
591 |
"The length returns the number of points in the LineString." |
|
|
592 |
return self.point_count |
|
|
593 |
|
|
|
594 |
@property |
|
|
595 |
def tuple(self): |
|
|
596 |
"Returns the tuple representation of this LineString." |
|
|
597 |
return tuple([self[i] for i in xrange(len(self))]) |
|
|
598 |
coords = tuple |
|
|
599 |
|
|
|
600 |
def _listarr(self, func): |
|
|
601 |
""" |
|
|
602 |
Internal routine that returns a sequence (list) corresponding with |
|
|
603 |
the given function. |
|
|
604 |
""" |
|
|
605 |
return [func(self.ptr, i) for i in xrange(len(self))] |
|
|
606 |
|
|
|
607 |
@property |
|
|
608 |
def x(self): |
|
|
609 |
"Returns the X coordinates in a list." |
|
|
610 |
return self._listarr(capi.getx) |
|
|
611 |
|
|
|
612 |
@property |
|
|
613 |
def y(self): |
|
|
614 |
"Returns the Y coordinates in a list." |
|
|
615 |
return self._listarr(capi.gety) |
|
29
|
616 |
|
|
0
|
617 |
@property |
|
|
618 |
def z(self): |
|
|
619 |
"Returns the Z coordinates in a list." |
|
|
620 |
if self.coord_dim == 3: |
|
|
621 |
return self._listarr(capi.getz) |
|
|
622 |
|
|
|
623 |
# LinearRings are used in Polygons. |
|
|
624 |
class LinearRing(LineString): pass |
|
|
625 |
|
|
|
626 |
class Polygon(OGRGeometry): |
|
|
627 |
|
|
|
628 |
def __len__(self): |
|
|
629 |
"The number of interior rings in this Polygon." |
|
|
630 |
return self.geom_count |
|
|
631 |
|
|
|
632 |
def __iter__(self): |
|
|
633 |
"Iterates through each ring in the Polygon." |
|
|
634 |
for i in xrange(self.geom_count): |
|
|
635 |
yield self[i] |
|
|
636 |
|
|
|
637 |
def __getitem__(self, index): |
|
|
638 |
"Gets the ring at the specified index." |
|
|
639 |
if index < 0 or index >= self.geom_count: |
|
|
640 |
raise OGRIndexError('index out of range: %s' % index) |
|
|
641 |
else: |
|
|
642 |
return OGRGeometry(capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs) |
|
|
643 |
|
|
|
644 |
# Polygon Properties |
|
|
645 |
@property |
|
|
646 |
def shell(self): |
|
|
647 |
"Returns the shell of this Polygon." |
|
|
648 |
return self[0] # First ring is the shell |
|
|
649 |
exterior_ring = shell |
|
|
650 |
|
|
|
651 |
@property |
|
|
652 |
def tuple(self): |
|
|
653 |
"Returns a tuple of LinearRing coordinate tuples." |
|
|
654 |
return tuple([self[i].tuple for i in xrange(self.geom_count)]) |
|
|
655 |
coords = tuple |
|
|
656 |
|
|
|
657 |
@property |
|
|
658 |
def point_count(self): |
|
|
659 |
"The number of Points in this Polygon." |
|
|
660 |
# Summing up the number of points in each ring of the Polygon. |
|
|
661 |
return sum([self[i].point_count for i in xrange(self.geom_count)]) |
|
|
662 |
|
|
|
663 |
@property |
|
|
664 |
def centroid(self): |
|
|
665 |
"Returns the centroid (a Point) of this Polygon." |
|
|
666 |
# The centroid is a Point, create a geometry for this. |
|
|
667 |
p = OGRGeometry(OGRGeomType('Point')) |
|
|
668 |
capi.get_centroid(self.ptr, p.ptr) |
|
|
669 |
return p |
|
|
670 |
|
|
|
671 |
# Geometry Collection base class. |
|
|
672 |
class GeometryCollection(OGRGeometry): |
|
|
673 |
"The Geometry Collection class." |
|
|
674 |
|
|
|
675 |
def __getitem__(self, index): |
|
|
676 |
"Gets the Geometry at the specified index." |
|
|
677 |
if index < 0 or index >= self.geom_count: |
|
|
678 |
raise OGRIndexError('index out of range: %s' % index) |
|
|
679 |
else: |
|
|
680 |
return OGRGeometry(capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs) |
|
29
|
681 |
|
|
0
|
682 |
def __iter__(self): |
|
|
683 |
"Iterates over each Geometry." |
|
|
684 |
for i in xrange(self.geom_count): |
|
|
685 |
yield self[i] |
|
|
686 |
|
|
|
687 |
def __len__(self): |
|
|
688 |
"The number of geometries in this Geometry Collection." |
|
|
689 |
return self.geom_count |
|
|
690 |
|
|
|
691 |
def add(self, geom): |
|
|
692 |
"Add the geometry to this Geometry Collection." |
|
|
693 |
if isinstance(geom, OGRGeometry): |
|
|
694 |
if isinstance(geom, self.__class__): |
|
|
695 |
for g in geom: capi.add_geom(self.ptr, g.ptr) |
|
|
696 |
else: |
|
|
697 |
capi.add_geom(self.ptr, geom.ptr) |
|
|
698 |
elif isinstance(geom, basestring): |
|
|
699 |
tmp = OGRGeometry(geom) |
|
|
700 |
capi.add_geom(self.ptr, tmp.ptr) |
|
|
701 |
else: |
|
|
702 |
raise OGRException('Must add an OGRGeometry.') |
|
|
703 |
|
|
|
704 |
@property |
|
|
705 |
def point_count(self): |
|
|
706 |
"The number of Points in this Geometry Collection." |
|
|
707 |
# Summing up the number of points in each geometry in this collection |
|
|
708 |
return sum([self[i].point_count for i in xrange(self.geom_count)]) |
|
|
709 |
|
|
|
710 |
@property |
|
|
711 |
def tuple(self): |
|
|
712 |
"Returns a tuple representation of this Geometry Collection." |
|
|
713 |
return tuple([self[i].tuple for i in xrange(self.geom_count)]) |
|
|
714 |
coords = tuple |
|
|
715 |
|
|
|
716 |
# Multiple Geometry types. |
|
|
717 |
class MultiPoint(GeometryCollection): pass |
|
|
718 |
class MultiLineString(GeometryCollection): pass |
|
|
719 |
class MultiPolygon(GeometryCollection): pass |
|
|
720 |
|
|
|
721 |
# Class mapping dictionary (using the OGRwkbGeometryType as the key) |
|
|
722 |
GEO_CLASSES = {1 : Point, |
|
|
723 |
2 : LineString, |
|
|
724 |
3 : Polygon, |
|
|
725 |
4 : MultiPoint, |
|
|
726 |
5 : MultiLineString, |
|
|
727 |
6 : MultiPolygon, |
|
|
728 |
7 : GeometryCollection, |
|
29
|
729 |
101: LinearRing, |
|
|
730 |
1 + OGRGeomType.wkb25bit : Point, |
|
|
731 |
2 + OGRGeomType.wkb25bit : LineString, |
|
|
732 |
3 + OGRGeomType.wkb25bit : Polygon, |
|
|
733 |
4 + OGRGeomType.wkb25bit : MultiPoint, |
|
|
734 |
5 + OGRGeomType.wkb25bit : MultiLineString, |
|
|
735 |
6 + OGRGeomType.wkb25bit : MultiPolygon, |
|
|
736 |
7 + OGRGeomType.wkb25bit : GeometryCollection, |
|
0
|
737 |
} |