1
0
mirror of https://github.com/django/django.git synced 2025-07-04 17:59:13 +00:00

gis: Added an intial ctypes interface to GDAL/OGR

(1) SpatialReference allows for rich access to properties of spatial 
      reference systems.
  (2) OGRGeometry provides access to the OGR Geometry classes -- may be
      accessed from models with the get_GEOM_ogr() extra instance method.


git-svn-id: http://code.djangoproject.com/svn/django/branches/gis@5397 bcc190cf-cafb-0310-a4f2-bffc1f526a37
This commit is contained in:
Justin Bronn 2007-06-01 01:15:35 +00:00
parent 38ff3cff45
commit b962a44a4d
27 changed files with 1674 additions and 31 deletions

View File

@ -1,5 +1,6 @@
# GEOS Routines # GEOS Routines
from django.contrib.gis.geos import GEOSGeometry, hex_to_wkt, centroid, area from django.contrib.gis.geos import GEOSGeometry, hex_to_wkt, centroid, area
from django.contrib.gis.gdal import OGRGeometry, SpatialReference
# Until model subclassing is a possibility, a mixin class is used to add # Until model subclassing is a possibility, a mixin class is used to add
# the necessary functions that may be contributed for geographic objects. # the necessary functions that may be contributed for geographic objects.
@ -12,6 +13,11 @@ class GeoMixin:
"Gets a GEOS Python object for the geometry." "Gets a GEOS Python object for the geometry."
return GEOSGeometry(getattr(self, field.attname), 'hex') return GEOSGeometry(getattr(self, field.attname), 'hex')
def _get_GEOM_ogr(self, field, srid):
"Gets an OGR Python object for the geometry."
return OGRGeometry(hex_to_wkt(getattr(self, field.attname)),
SpatialReference('EPSG:%d' % srid))
def _get_GEOM_wkt(self, field): def _get_GEOM_wkt(self, field):
"Gets the WKT of the geometry." "Gets the WKT of the geometry."
hex = getattr(self, field.attname) hex = getattr(self, field.attname)

View File

@ -6,7 +6,7 @@ from django.utils.functional import curry
from django.contrib.gis.geos import GEOSGeometry, GEOSException from django.contrib.gis.geos import GEOSGeometry, GEOSException
#TODO: Flesh out widgets. #TODO: Flesh out widgets.
#TODO: geos operations through fields as proxy. #TODO: GEOS and GDAL/OGR operations through fields as proxy.
#TODO: pythonic usage, like "for point in zip.polygon" and "if point in polygon". #TODO: pythonic usage, like "for point in zip.polygon" and "if point in polygon".
class GeometryField(Field): class GeometryField(Field):
@ -38,8 +38,7 @@ class GeometryField(Field):
AddGeometryColumn(...) PostGIS (and OGC standard) function. AddGeometryColumn(...) PostGIS (and OGC standard) function.
Takes the style object (provides syntax highlighting) as well as the Takes the style object (provides syntax highlighting) as well as the
database table and field. The dimensions can be specified via database table and field.
the dim keyword as well.
""" """
sql = style.SQL_KEYWORD('SELECT ') + \ sql = style.SQL_KEYWORD('SELECT ') + \
style.SQL_TABLE('AddGeometryColumn') + '(' + \ style.SQL_TABLE('AddGeometryColumn') + '(' + \
@ -81,8 +80,9 @@ class GeometryField(Field):
def contribute_to_class(self, cls, name): def contribute_to_class(self, cls, name):
super(GeometryField, self).contribute_to_class(cls, name) super(GeometryField, self).contribute_to_class(cls, name)
# Adding the WKT accessor function for geometry # Adding needed accessor functions
setattr(cls, 'get_%s_geos' % self.name, curry(cls._get_GEOM_geos, field=self)) setattr(cls, 'get_%s_geos' % self.name, curry(cls._get_GEOM_geos, field=self))
setattr(cls, 'get_%s_ogr' % self.name, curry(cls._get_GEOM_ogr, field=self, srid=self._srid))
setattr(cls, 'get_%s_wkt' % self.name, curry(cls._get_GEOM_wkt, field=self)) setattr(cls, 'get_%s_wkt' % self.name, curry(cls._get_GEOM_wkt, field=self))
setattr(cls, 'get_%s_centroid' % self.name, curry(cls._get_GEOM_centroid, field=self)) setattr(cls, 'get_%s_centroid' % self.name, curry(cls._get_GEOM_centroid, field=self))
setattr(cls, 'get_%s_area' % self.name, curry(cls._get_GEOM_area, field=self)) setattr(cls, 'get_%s_area' % self.name, curry(cls._get_GEOM_area, field=self))

View File

@ -1,6 +1,5 @@
# This module is meant to re-define the helper routines used by the # This module is meant to re-define the helper routines used by the
# django.db.models.query objects to be customized for PostGIS. # django.db.models.query objects to be customized for PostGIS.
from copy import copy
from django.db import backend from django.db import backend
from django.db.models.query import LOOKUP_SEPARATOR, find_field, FieldFound, QUERY_TERMS, get_where_clause from django.db.models.query import LOOKUP_SEPARATOR, find_field, FieldFound, QUERY_TERMS, get_where_clause
from django.utils.datastructures import SortedDict from django.utils.datastructures import SortedDict

View File

@ -0,0 +1,103 @@
# types and ctypes
from types import StringType
from ctypes import c_char_p, c_int, c_void_p, byref, string_at
# The GDAL C library, OGR exceptions, and the Layer object.
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.OGRError import OGRException, check_err
from django.contrib.gis.gdal.Layer import Layer
"""
DataSource is a wrapper for the OGR Data Source object, which provides
an interface for reading vector geometry data from many different file
formats (including ESRI shapefiles).
Example:
ds = DataSource('/home/foo/bar.shp')
for layer in ds:
for feature in layer:
# Getting the geometry for the feature.
g = feature.geom
# Getting the 'description' field for the feature.
desc = feature['description']
More documentation forthcoming.
"""
# For more information, see the OGR C API source code:
# http://www.gdal.org/ogr/ogr__api_8h.html
#
# The OGR_DS* routines are relevant here.
class DataSource(object):
"Wraps an OGR Data Source object."
_ds = 0 # Initially NULL
#### Python 'magic' routines ####
def __init__(self, ds_file):
# Registering all the drivers, this needs to be done
# _before_ we try to open up a data source.
if not lgdal.OGRRegisterAll():
raise OGRException, 'Could not register all data source drivers!'
# The data source driver is a void pointer.
ds_driver = c_void_p()
# OGROpen will auto-detect the data source type.
ds = lgdal.OGROpen(c_char_p(ds_file), c_int(0), byref(ds_driver))
# Raise an exception if the returned pointer is NULL
if not ds:
self._ds = False
raise OGRException, 'Invalid data source file "%s"' % ds_file
else:
self._ds = ds
self._driver = ds_driver
def __del__(self):
"This releases the reference to the data source (destroying it if it's the only one)."
if self._ds: lgdal.OGRReleaseDataSource(self._ds)
def __iter__(self):
"Allows for iteration over the layers in a data source."
for i in xrange(self.layer_count):
yield self.__getitem__(i)
def __getitem__(self, index):
"Allows use of the index [] operator to get a layer at the index."
if isinstance(index, StringType):
l = lgdal.OGR_DS_GetLayerByName(self._ds, c_char_p(index))
if not l: raise IndexError, 'invalid OGR Layer name given: "%s"' % index
else:
if index < 0 or index >= self.layer_count:
raise IndexError, 'index out of range'
l = lgdal.OGR_DS_GetLayer(self._ds, c_int(index))
return Layer(l)
def __len__(self):
"Returns the number of layers within the data source."
return self.layer_count
def __str__(self):
"Returns OGR GetName and Driver for the Data Source."
return '%s (%s)' % (self.name, self.driver)
#### DataSource Properties ####
@property
def driver(self):
"Returns the name of the data source driver."
return string_at(lgdal.OGR_Dr_GetName(self._driver))
@property
def layer_count(self):
"Returns the number of layers in the data source."
return lgdal.OGR_DS_GetLayerCount(self._ds)
@property
def name(self):
"Returns the name of the data source."
return string_at(lgdal.OGR_DS_GetName(self._ds))

View File

@ -0,0 +1,104 @@
# types and ctypes
import types
from ctypes import c_char_p, c_int, string_at
# The GDAL C library, OGR exception, and the Field object
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.OGRError import OGRException
from django.contrib.gis.gdal.Field import Field
from django.contrib.gis.gdal.OGRGeometry import OGRGeometry, OGRGeomType
# For more information, see the OGR C API source code:
# http://www.gdal.org/ogr/ogr__api_8h.html
#
# The OGR_F_* routines are relevant here.
class Feature(object):
"A class that wraps an OGR Feature, needs to be instantiated from a Layer object."
_feat = 0 # Initially NULL
#### Python 'magic' routines ####
def __init__(self, f):
"Needs a C pointer (Python integer in ctypes) in order to initialize."
if not f:
raise OGRException, 'Cannot create OGR Feature, invalid pointer given.'
self._feat = f
self._fdefn = lgdal.OGR_F_GetDefnRef(f)
def __del__(self):
"Releases a reference to this object."
if self._fdefn: lgdal.OGR_FD_Release(self._fdefn)
def __getitem__(self, index):
"Gets the Field at the specified index."
if isinstance(index, types.StringType):
i = self.index(index)
else:
if index < 0 or index > self.num_fields:
raise IndexError, 'index out of range'
i = index
return Field(lgdal.OGR_F_GetFieldDefnRef(self._feat, c_int(i)),
string_at(lgdal.OGR_F_GetFieldAsString(self._feat, c_int(i))))
def __iter__(self):
"Iterates over each field in the Feature."
for i in xrange(self.num_fields):
yield self.__getitem__(i)
def __len__(self):
"Returns the count of fields in this feature."
return self.num_fields
def __str__(self):
"The string name of the feature."
return 'Feature FID %d in Layer<%s>' % (self.fid, self.layer_name)
def __eq__(self, other):
"Does equivalence testing on the features."
if lgdal.OGR_F_Equal(self._feat, other._feat):
return True
else:
return False
#### Feature Properties ####
@property
def fid(self):
"Returns the feature identifier."
return lgdal.OGR_F_GetFID(self._feat)
@property
def layer_name(self):
"Returns the name of the layer for the feature."
return string_at(lgdal.OGR_FD_GetName(self._fdefn))
@property
def num_fields(self):
"Returns the number of fields in the Feature."
return lgdal.OGR_F_GetFieldCount(self._feat)
@property
def fields(self):
"Returns a list of fields in the Feature."
return [ string_at(lgdal.OGR_Fld_GetNameRef(lgdal.OGR_FD_GetFieldDefn(self._fdefn, i)))
for i in xrange(self.num_fields) ]
@property
def geom(self):
"Returns the OGR Geometry for this Feature."
# A clone is used, so destruction of the Geometry won't bork the Feature.
return OGRGeometry(lgdal.OGR_G_Clone(lgdal.OGR_F_GetGeometryRef(self._feat)))
@property
def geom_type(self):
"Returns the OGR Geometry Type for this Feture."
return OGRGeomType(lgdal.OGR_FD_GetGeomType(self._fdefn))
#### Feature Methods ####
def index(self, field_name):
"Returns the index of the given field name."
i = lgdal.OGR_F_GetFieldIndex(self._feat, c_char_p(field_name))
if i < 0: raise IndexError, 'invalid OFT field name given: "%s"' % field_name
return i
def clone(self):
"Clones this Feature."
return Feature(lgdal.OGR_F_Clone(self._feat))

View File

@ -0,0 +1,83 @@
from ctypes import string_at
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.OGRError import OGRException
# For more information, see the OGR C API source code:
# http://www.gdal.org/ogr/ogr__api_8h.html
#
# The OGR_Fld_* routines are relevant here.
class Field(object):
"A class that wraps an OGR Field."
_fld = 0 # Initially NULL
#### Python 'magic' routines ####
def __init__(self, fld, val=''):
"Needs a C pointer (Python integer in ctypes) in order to initialize."
if not fld:
raise OGRException, 'Cannot create OGR Field, invalid pointer given.'
self._fld = fld
self._val = val
# Setting the class depending upon the OGR Field Type (OFT)
self.__class__ = FIELD_CLASSES[self.type]
def __str__(self):
"Returns the string representation of the Field."
return '%s (%s)' % (self.name, self.__class__.__name__)
#### Field Properties ####
@property
def name(self):
"Returns the name of the field."
return string_at(lgdal.OGR_Fld_GetNameRef(self._fld))
@property
def type(self):
"Returns the type of this field."
return lgdal.OGR_Fld_GetType(self._fld)
@property
def value(self):
"Returns the value of this type of field."
return self._val
# The Field sub-classes for each OGR Field type.
class OFTInteger(Field):
@property
def value(self):
"Returns an integer contained in this field."
return int(self._val)
class OFTIntegerList(Field): pass
class OFTReal(Field):
@property
def value(self):
"Returns a float contained in this field."
return float(self._val)
class OFTRealList(Field): pass
class OFTString(Field): pass
class OFTStringList(Field): pass
class OFTWideString(Field): pass
class OFTWideStringList(Field): pass
class OFTBinary(Field): pass
class OFTDate(Field): pass
class OFTTime(Field): pass
class OFTDateTime(Field): pass
# Class mapping dictionary for OFT Types
FIELD_CLASSES = { 0 : OFTInteger,
1 : OFTIntegerList,
2 : OFTReal,
3 : OFTRealList,
4 : OFTString,
5 : OFTStringList,
6 : OFTWideString,
7 : OFTWideStringList,
8 : OFTBinary,
9 : OFTDate,
10 : OFTTime,
11 : OFTDateTime,
}

View File

@ -0,0 +1,28 @@
Copyright (c) 2007, Justin Bronn
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of OGRGeometry nor the names of its contributors may be used
to endorse or promote products derived from this software without
specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

View File

@ -0,0 +1,90 @@
# Needed ctypes routines
from ctypes import c_int, c_long, c_void_p, string_at
# The GDAL C Library
from django.contrib.gis.gdal.libgdal import lgdal
# Other GDAL imports.
from django.contrib.gis.gdal.Feature import Feature
from django.contrib.gis.gdal.OGRError import OGRException
from django.contrib.gis.gdal.SpatialReference import SpatialReference
# For more information, see the OGR C API source code:
# http://www.gdal.org/ogr/ogr__api_8h.html
#
# The OGR_L_* routines are relevant here.
get_srs = lgdal.OGR_L_GetSpatialRef
get_srs.restype = c_void_p
get_srs.argtypes = [c_void_p]
class Layer(object):
"A class that wraps an OGR Layer, needs to be instantiated from a DataSource object."
_layer = 0 # Initially NULL
#### Python 'magic' routines ####
def __init__(self, l):
"Needs a C pointer (Python/ctypes integer) in order to initialize."
if not l:
raise OGRException, 'Cannot create Layer, invalid pointer given'
self._layer = l
self._ldefn = lgdal.OGR_L_GetLayerDefn(l)
def __getitem__(self, index):
"Gets the Feature at the specified index."
if index < 0 or index >= self.num_feat:
raise IndexError, 'index out of range'
return Feature(lgdal.OGR_L_GetFeature(self._layer, c_long(index)))
def __iter__(self):
"Iterates over each Feature in the Layer."
# Resetting the Layer before beginning iteration
lgdal.OGR_L_ResetReading(self._layer)
# Incrementing over each feature in the layer, and yielding
# to the caller of the function.
for i in xrange(self.num_feat):
yield self.__getitem__(i)
def __len__(self):
"The length is the number of features."
return self.num_feat
def __str__(self):
"The string name of the layer."
return self.name
#### Layer properties ####
@property
def name(self):
"Returns the name of this layer in the Data Source."
return string_at(lgdal.OGR_FD_GetName(self._ldefn))
@property
def num_feat(self, force=1):
"Returns the number of features in the Layer."
return lgdal.OGR_L_GetFeatureCount(self._layer, c_int(force))
@property
def num_fields(self):
"Returns the number of fields in the Layer."
return lgdal.OGR_FD_GetFieldCount(self._ldefn)
@property
def geom_type(self):
"Returns the geometry type (OGRGeomType) of the Layer."
return lgdal.OGR_FD_GetGeomType(self._ldefn)
@property
def srs(self):
"Returns the Spatial Reference used in this Layer."
return SpatialReference(lgdal.OSRClone(lgdal.OGR_L_GetSpatialRef(self._layer)), 'ogr')
@property
def fields(self):
"Returns a list of the fields available in this Layer."
return [ string_at(lgdal.OGR_Fld_GetNameRef(lgdal.OGR_FD_GetFieldDefn(self._ldefn, i)))
for i in xrange(self.num_fields) ]

View File

@ -0,0 +1,35 @@
# OGR Error Codes
OGRERR_NONE = 0
OGRERR_NOT_ENOUGH_DATA = 1
OGRERR_NOT_ENOUGH_MEMORY = 2
OGRERR_UNSUPPORTED_GEOMETRY_TYPE = 3
OGRERR_UNSUPPORTED_OPERATION = 4
OGRERR_CORRUPT_DATA = 5
OGRERR_FAILURE = 6
OGRERR_UNSUPPORTED_SRS = 7
# OGR & SRS Exceptions
class OGRException(Exception): pass
class SRSException(Exception): pass
def check_err(code, msg=False):
"Checks the given OGRERR, and raises an exception where appropriate."
if code == OGRERR_NONE:
return
elif code == OGRERR_NOT_ENOUGH_DATA:
raise OGRException, 'Not enough data!'
elif code == OGRERR_NOT_ENOUGH_MEMORY:
raise OGRException, 'Not enough memory!'
elif code == OGRERR_UNSUPPORTED_GEOMETRY_TYPE:
raise OGRException, 'Unsupported Geometry Type!'
elif code == OGRERR_UNSUPPORTED_OPERATION:
raise OGRException, 'Unsupported Operation!'
elif code == OGRERR_CORRUPT_DATA:
raise OGRException, 'Corrupt Data!'
elif code == OGRERR_FAILURE:
raise OGRException, 'OGR Failure!'
elif code == OGRERR_UNSUPPORTED_SRS:
raise SRSException, 'Unsupported SRS!'
else:
raise OGRException, 'Unknown error code: "%s"' % str(code)

View File

@ -0,0 +1,340 @@
# types & ctypes
from types import StringType
from ctypes import \
byref, string_at, create_string_buffer, POINTER, \
c_char_p, c_double, c_int, c_void_p
# Getting the GDAL C library and error checking facilities
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.OGRError import check_err, OGRException
from django.contrib.gis.gdal.SpatialReference import SpatialReference, CoordTransform
# 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.
#### ctypes prototypes ####
def pnt_func(f):
"For accessing point information."
f.restype = c_double
f.argtypes = [c_void_p, c_int]
return f
getx = pnt_func(lgdal.OGR_G_GetX)
gety = pnt_func(lgdal.OGR_G_GetY)
getz = pnt_func(lgdal.OGR_G_GetZ)
#### OGRGeomType ####
class OGRGeomType(object):
"Encapulates OGR Geometry Types."
# Ordered array of acceptable strings and their corresponding OGRwkbGeometryType
__ogr_str = ['Unknown', 'Point', 'LineString', 'Polygon', 'MultiPoint',
'MultiLineString', 'MultiPolygon', 'GeometryCollection',
'None', 'LinearRing']
__ogr_int = [0, 1, 2, 3, 4, 5, 6, 7, 100, 101]
def __init__(self, input):
"Figures out the correct OGR Type based upon the input."
if isinstance(input, OGRGeomType):
self._index = input._index
elif isinstance(input, StringType):
idx = self._has_str(self.__ogr_str, input)
if not idx:
raise OGRException, 'Invalid OGR String Type "%s"' % input
self._index = idx
elif isinstance(input, int):
if not input in self.__ogr_int:
raise OGRException, 'Invalid OGR Integer Type: %d' % input
self._index = self.__ogr_int.index(input)
else:
raise TypeError, 'Invalid OGR Input type given!'
def __str__(self):
"Returns a short-hand string form of the OGR Geometry type."
return self.__ogr_str[self._index]
def __eq__(self, other):
"""Does an equivalence test on the OGR type with the given
other OGRGeomType, the short-hand string, or the integer."""
if isinstance(other, OGRGeomType):
return self._index == other._index
elif isinstance(other, StringType):
idx = self._has_str(self.__ogr_str, other)
if idx: return self._index == idx
return False
elif isinstance(other, int):
if not other in self.__ogr_int: return False
return self.__ogr_int.index(other) == self._index
else:
raise TypeError, 'Cannot compare with type: %s' % str(type(other))
def _has_str(self, arr, s):
slow = s.lower()
for i in xrange(len(arr)):
if slow == arr[i].lower(): return i
return None
@property
def django(self):
"Returns the Django GeometryField for this OGR Type."
s = self.__ogr_str[self._index]
if s in ('Unknown', 'LinearRing'):
return None
else:
return s + 'Field'
@property
def num(self):
"Returns the OGRwkbGeometryType number for the OGR Type."
return self.__ogr_int[self._index]
#### OGRGeometry Class ####
class OGRGeometry(object):
"Generally encapsulates an OGR geometry."
_g = 0 # Initially NULL
def __init__(self, input, srs=False):
"Initializes Geometry on either WKT or an OGR pointer as input."
if isinstance(input, StringType):
# First, trying the input as WKT
buf = c_char_p(input)
g = c_void_p()
# Getting the spatial
if not isinstance(srs, SpatialReference):
s = SpatialReference() # creating an empty spatial reference
else:
s = srs.clone() # cloning the given spatial reference
try:
check_err(lgdal.OGR_G_CreateFromWkt(byref(buf), s._srs, byref(g)))
except OGRException, msg:
try:
ogr_t = OGRGeomType(input) # Seeing if the input is a valid short-hand string
g = lgdal.OGR_G_CreateGeometry(ogr_t.num)
except:
raise OGRException, 'Could not initialize on WKT "%s"' % input
elif isinstance(input, OGRGeomType):
g = lgdal.OGR_G_CreateGeometry(input.num)
elif isinstance(input, int):
# OGR Pointer (integer) was the input
g = input
else:
raise OGRException, 'Type of input cannot be determined!'
# Now checking the Geometry pointer before finishing initialization
if not g:
raise OGRException, 'Cannot create OGR Geometry from input: %s' % str(input)
self._g = g
# 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._g: lgdal.OGR_G_DestroyGeometry(self._g)
def __eq__(self, other):
"Is this Geometry equal to the other?"
return lgdal.OGR_G_Equals(self._g, other._g)
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 lgdal.OGR_G_GetDimension(self._g)
@property
def coord_dim(self):
"Returns the coordinate dimension of the Geometry."
return lgdal.OGR_G_GetCoordinateDimension(self._g)
@property
def geom_count(self):
"The number of elements in this Geometry."
return lgdal.OGR_G_GetGeometryCount(self._g)
@property
def point_count(self):
"The number of Points in this Geometry."
return lgdal.OGR_G_GetPointCount(self._g)
@property
def srs(self):
"Returns the Spatial Reference for this Geometry."
return SpatialReference(lgdal.OSRClone(lgdal.OGR_G_GetSpatialReference(self._g)), 'ogr')
@property
def geom_type(self):
"Returns the Type for this Geometry."
return OGRGeomType(lgdal.OGR_G_GetGeometryType(self._g))
@property
def geom_name(self):
"Returns the Name of this Geometry."
return string_at(lgdal.OGR_G_GetGeometryName(self._g))
@property
def wkt(self):
"Returns the WKT form of the Geometry."
buf = c_char_p()
check_err(lgdal.OGR_G_ExportToWkt(self._g, byref(buf)))
return string_at(buf)
#### Geometry Methods ####
def clone(self):
"Clones this OGR Geometry."
return OGRGeometry(lgdal.OGR_G_Clone(self._g))
def transform(self, coord_trans):
"Transforms this Geometry with the given CoordTransform object."
if not isinstance(coord_trans, CoordTransform):
raise OGRException, 'CoordTransform object required for transform.'
check_err(lgdal.OGR_G_Transform(self._g, coord_trans._ct))
def transform_to(self, srs):
"Transforms this Geometry with the given SpatialReference."
if not isinstance(srs, SpatialReference):
raise OGRException, 'SpatialReference object required for transform_to.'
check_err(lgdal.OGR_G_TransformTo(self._g, srs._srs))
# The subclasses for OGR Geometry.
class Point(OGRGeometry):
@property
def x(self):
"Returns the X coordinate for this Point."
return getx(self._g, c_int(0))
@property
def y(self):
"Returns the Y coordinate for this Point."
return gety(self._g, c_int(0))
@property
def z(self):
"Returns the Z coordinate for this Point."
return getz(self._g, c_int(0))
@property
def tuple(self):
"Returns the tuple of this point."
if self.coord_dim == 1:
return (self.x,)
elif self.coord_dim == 2:
return (self.x, self.y)
elif self.coord_dim == 3:
return (self.x, self.y, self.z)
class LineString(OGRGeometry):
def __getitem__(self, index):
"Returns the Point at the given index."
if index > 0 or index < self.point_count:
x = c_double()
y = c_double()
z = c_double()
lgdal.OGR_G_GetPoint(self._g, c_int(index),
byref(x), byref(y), byref(z))
if self.coord_dim == 1:
return (x.value,)
elif self.coord_dim == 2:
return (x.value, y.value)
elif self.coord_dim == 3:
return (x.value, y.value, z.value)
else:
raise IndexError, 'index out of range'
def __iter__(self):
"Iterates over each point in the LineString."
for i in xrange(self.point_count):
yield self.__getitem__(i)
def __len__(self, index):
"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.__getitem__(i) for i in xrange(self.point_count))
# 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.__getitem__(i)
def __getitem__(self, index):
"Gets the ring at the specified index."
if index < 0 or index >= self.geom_count:
raise IndexError, 'index out of range'
else:
return OGRGeometry(lgdal.OGR_G_Clone(lgdal.OGR_G_GetGeometryRef(self._g, c_int(index))))
# Polygon Properties
@property
def shell(self):
"Returns the shell of this Polygon."
return self.__getitem__(0) # First ring is the shell
@property
def tuple(self):
"Returns a tuple of LinearRing coordinate tuples."
return tuple(self.__getitem__(i).tuple for i in xrange(self.geom_count))
# 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 IndexError, 'index out of range'
else:
return OGRGeometry(lgdal.OGR_G_Clone(lgdal.OGR_G_GetGeometryRef(self._g, c_int(index))))
def __iter__(self):
"Iterates over each Geometry."
for i in xrange(self.geom_count):
yield self.__getitem__(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 not isinstance(geom, OGRGeometry):
raise OGRException, 'Must add an OGRGeometry.'
lgdal.OGR_G_AddGeometry(self._g, geom._g)
# 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,
}

View File

@ -0,0 +1,299 @@
# Getting what we need from ctypes
import re
from types import StringType, TupleType
from ctypes import \
c_char_p, c_int, c_double, c_void_p, POINTER, \
byref, string_at, create_string_buffer
# Getting the GDAL C Library
from django.contrib.gis.gdal.libgdal import lgdal
# Getting the error checking routine and exceptions
from django.contrib.gis.gdal.OGRError import check_err, OGRException, SRSException
#### ctypes function prototypes ####
def ellipsis_func(f):
"""Creates a ctypes function prototype for OSR ellipsis property functions,
e.g., OSRGetSemiMajor, OSRGetSemiMinor, OSRGetInvFlattening."""
f.restype = c_double
f.argtypes = [c_void_p, POINTER(c_int)]
return f
# Getting the semi_major, semi_minor, and flattening functions.
semi_major = ellipsis_func(lgdal.OSRGetSemiMajor)
semi_minor = ellipsis_func(lgdal.OSRGetSemiMinor)
invflattening = ellipsis_func(lgdal.OSRGetInvFlattening)
def units_func(f):
"""Creates a ctypes function prototype for OSR units functions,
e.g., OSRGetAngularUnits, OSRGetLinearUnits."""
f.restype = c_double
f.argtypes = [c_void_p, POINTER(c_char_p)]
return f
# Getting the angular_units, linear_units functions
linear_units = units_func(lgdal.OSRGetLinearUnits)
angular_units = units_func(lgdal.OSRGetAngularUnits)
#### Spatial Reference class. ####
class SpatialReference(object):
"""A wrapper for the OGRSpatialReference object. According to the GDAL website,
the SpatialReference object 'provide[s] services to represent coordinate systems
(projections and datums) and to transform between them.'"""
_srs = 0 # Initially NULL
# Well-Known Geographical Coordinate System Name
_well_known = {'WGS84':4326, 'WGS72':4322, 'NAD27':4267, 'NAD83':4269}
_epsg_regex = re.compile('^EPSG:(?P<epsg>\d+)$', re.I)
#### Python 'magic' routines ####
def __init__(self, input='', srs_type='wkt'):
"Creates a spatial reference object from the given OGC Well Known Text (WKT)."
# Creating an initial empty string buffer.
buf = c_char_p('')
if isinstance(input, StringType):
# Is this an EPSG well known name?
m = self._epsg_regex.match(input)
if m:
srs_type = 'epsg'
input = int(m.group('epsg'))
# Is this a short-hand well known name?
elif input in self._well_known:
srs_type = 'epsg'
input = self._well_known[input]
elif srs_type == 'proj':
pass
else:
buf = c_char_p(input)
elif isinstance(input, int):
if srs_type not in ('epsg', 'ogr'):
raise SRSException, 'Integer input requires SRS type of "ogr" or "epsg".'
else:
raise TypeError, 'Invalid SRS type "%s"' % srs_type
# Calling OSRNewSpatialReference with the string buffer.
if srs_type == 'ogr':
srs = input # Input is OGR pointer
else:
srs = lgdal.OSRNewSpatialReference(buf)
# If the pointer is NULL, throw an exception.
if not srs:
raise SRSException, 'Could not create spatial reference from WKT!'
else:
self._srs = srs
# Post-processing if in PROJ.4 or EPSG formats.
if srs_type == 'proj': self.import_proj(input)
elif srs_type == 'epsg': self.import_epsg(input)
def __del__(self):
"Destroys this spatial reference."
if self._srs: lgdal.OSRRelease(self._srs)
def __getitem__(self, target):
"""Returns the value of the given string attribute node, None if the node doesn't exist.
Can also take a tuple as a parameter, (target, child), where child is the child index to get."""
if isinstance(target, TupleType):
return self.attr_value(*target)
else:
return self.attr_value(target)
def __str__(self):
"The string representation uses 'pretty' WKT."
return self.pretty_wkt
def _string_ptr(self, ptr):
"Returns the string at the pointer if it is valid, None if the pointer is NULL."
if not ptr: return None
else: return string_at(ptr)
#### SpatialReference Methods ####
def auth_name(self, target):
"Getting the authority name for the target node."
ptr = lgdal.OSRGetAuthorityName(self._srs, c_char_p(target))
return self._string_ptr(ptr)
def auth_code(self, target):
"Getting the authority code for the given target node."
ptr = lgdal.OSRGetAuthorityCode(self._srs, c_char_p(target))
return self._string_ptr(ptr)
def attr_value(self, target, index=0):
"""The attribute value for the given target node (e.g. 'PROJCS'). The index keyword
specifies an index of the child node to return."""
ptr = lgdal.OSRGetAttrValue(self._srs, c_char_p(target), c_int(index))
return self._string_ptr(ptr)
def validate(self):
"Checks to see if the given spatial reference is valid."
check_err(lgdal.OSRValidate(self._srs))
def clone(self):
"Returns a clone of this Spatial Reference."
return SpatialReference(lgdal.OSRClone(self._srs), 'ogr')
@property
def name(self):
"Returns the name of this Spatial Reference."
if self.projected: return self.attr_value('PROJCS')
elif self.geographic: return self.attr_value('GEOGCS')
elif self.local: return self.attr_value('LOCAL_CS')
else: return None
#### Unit Properties ####
def _cache_linear(self):
"Caches the linear units value and name."
if not hasattr(self, '_linear_units') or not hasattr(self, '_linear_name'):
name_buf = c_char_p()
self._linear_units = linear_units(self._srs, byref(name_buf))
self._linear_name = string_at(name_buf)
@property
def linear_name(self):
"Returns the name of the linear units."
self._cache_linear()
return self._linear_name
@property
def linear_units(self):
"Returns the value of the linear units."
self._cache_linear()
return self._linear_units
def _cache_angular(self):
"Caches the angular units value and name."
name_buf = c_char_p()
if not hasattr(self, '_angular_units') or not hasattr(self, '_angular_name'):
self._angular_units = angular_units(self._srs, byref(name_buf))
self._angular_name = string_at(name_buf)
@property
def angular_name(self):
"Returns the name of the angular units."
self._cache_angular()
return self._angular_name
@property
def angular_units(self):
"Returns the value of the angular units."
self._cache_angular()
return self._angular_units
#### Spheroid/Ellipsis Properties ####
@property
def semi_major(self):
"Gets the Semi Major Axis for this Spatial Reference."
err = c_int(0)
sm = semi_major(self._srs, byref(err))
check_err(err.value)
return sm
@property
def semi_minor(self):
"Gets the Semi Minor Axis for this Spatial Reference."
err = c_int()
sm = semi_minor(self._srs, byref(err))
check_err(err.value)
return sm
@property
def inverse_flattening(self):
"Gets the Inverse Flattening for this Spatial Reference."
err = c_int()
inv_flat = invflattening(self._srs, byref(err))
check_err(err.value)
return inv_flat
#### Boolean Properties ####
@property
def geographic(self):
"Returns True if this SpatialReference is geographic (root node is GEOGCS)."
if lgdal.OSRIsGeographic(self._srs): return True
else: return False
@property
def local(self):
"Returns True if this SpatialReference is local (root node is LOCAL_CS)."
if lgdal.OSRIsLocal(self._srs): return True
else: return False
@property
def projected(self):
"Returns True if this SpatialReference is a projected coordinate system (root node is PROJCS)."
if lgdal.OSRIsProjected(self._srs): return True
else: return False
#### Import Routines #####
def import_wkt(self, wkt):
"Imports the Spatial Reference from OGC WKT (string)"
buf = create_string_buffer(wkt)
check_err(lgdal.OSRImportFromWkt(self._srs, byref(buf)))
def import_proj(self, proj):
"Imports the Spatial Reference from a PROJ.4 string."
check_err(lgdal.OSRImportFromProj4(self._srs, create_string_buffer(proj)))
def import_epsg(self, epsg):
"Imports the Spatial Reference from the EPSG code (an integer)."
check_err(lgdal.OSRImportFromEPSG(self._srs, c_int(epsg)))
def import_xml(self, xml):
"Imports the Spatial Reference from an XML string."
check_err(lgdal.OSRImportFromXML(self._srs, create_string_buffer(xml)))
#### Export Properties ####
@property
def wkt(self):
"Returns the WKT representation of this Spatial Reference."
w = c_char_p()
check_err(lgdal.OSRExportToWkt(self._srs, byref(w)))
return string_at(w)
@property
def pretty_wkt(self, simplify=0):
"Returns the 'pretty' representation of the WKT."
w = c_char_p()
check_err(lgdal.OSRExportToPrettyWkt(self._srs, byref(w), c_int(simplify)))
return string_at(w)
@property
def proj(self):
"Returns the PROJ.4 representation for this Spatial Reference."
w = c_char_p()
check_err(lgdal.OSRExportToProj4(self._srs, byref(w)))
return string_at(w)
@property
def xml(self, dialect=''):
"Returns the XML representation of this Spatial Reference."
w = c_char_p()
check_err(lgdal.OSRExportToXML(self._srs, byref(w), create_string_buffer(dialect)))
return string_at(w)
class CoordTransform(object):
"A coordinate system transformation object."
_ct = 0 # Initially NULL
def __init__(self, source, target):
"Initializes on a source and target SpatialReference objects."
if not isinstance(source, SpatialReference) or not isinstance(target, SpatialReference):
raise SRSException, 'source and target must be of type SpatialReference'
ct = lgdal.OCTNewCoordinateTransformation(source._srs, target._srs)
if not ct:
raise SRSException, 'could not intialize CoordTransform object'
self._ct = ct
self._srs1_name = source.name
self._srs2_name = target.name
def __del__(self):
"Deletes this Coordinate Transformation object."
if self._ct: lgdal.OCTDestroyCoordinateTransformation(self._ct)
def __str__(self):
return 'Transform from "%s" to "%s"' % (str(self._srs1_name), str(self._srs2_name))

View File

@ -0,0 +1,5 @@
from DataSource import DataSource
from SpatialReference import SpatialReference, CoordTransform
from OGRGeometry import OGRGeometry, OGRGeomType
from OGRError import check_err, OGRException, SRSException

View File

@ -0,0 +1,22 @@
import os, sys
from ctypes import CDLL
if os.name == 'nt':
# Windows NT library
lib_name = 'libgdal-1.dll'
elif os.name == 'posix':
platform = os.uname()[0]
if platform == 'Linux':
# Linux shared library
lib_name = 'libgdal.so'
elif platform == 'Darwin':
# Mac OSX Shared Library
lib_name = 'libgdal.dylib'
else:
raise GDALException, 'Unknown POSIX platform "%s"' % platform
else:
raise GDALException, 'Unsupported OS "%s"' % os.name
# The GDAL C library
lgdal = CDLL(lib_name)

View File

@ -1,26 +1,29 @@
import re import re
from django.db import models from django.db import models
# Checking for the presence of GDAL
try:
from django.contrib.gis.gdal import SpatialReference
HAS_OSR = True
except ImportError:
HAS_OSR = False
""" """
Models for the PostGIS/OGC database tables. Models for the PostGIS/OGC database tables.
""" """
# For pulling out the spheroid from the spatial reference string. # For pulling out the spheroid from the spatial reference string. This
# regular expression is used only if the user does not have GDAL installed.
# TODO: Flattening not used in all ellipsoids, could also be a minor axis, or 'b' # TODO: Flattening not used in all ellipsoids, could also be a minor axis, or 'b'
# parameter. # parameter.
spheroid_regex = re.compile(r'.+SPHEROID\[\"(?P<name>.+)\",(?P<major>\d+(\.\d+)?),(?P<flattening>\d{3}\.\d+),') spheroid_regex = re.compile(r'.+SPHEROID\[\"(?P<name>.+)\",(?P<major>\d+(\.\d+)?),(?P<flattening>\d{3}\.\d+),')
# For pulling out the projected coordinate system units. Python regexs are greedy
# by default, so this should get the units of the projection instead (PROJCS)
# of the units for the geographic coordinate system (GEOGCS).
unit_regex = re.compile(r'^PROJCS.+UNIT\[\"(?P<units>[a-z]+)\", ?(?P<conversion>[0-9\.]+), ?(AUTHORITY\[\"(?P<authority>[a-z0-9 \.]+)\",[ ]?\"(?P<code>\d+)\"\])?', re.I)
# This is the global 'geometry_columns' from PostGIS. # This is the global 'geometry_columns' from PostGIS.
# See PostGIS Documentation at Ch. 4.2.2 # See PostGIS Documentation at Ch. 4.2.2
class GeometryColumns(models.Model): class GeometryColumns(models.Model):
f_table_catalog = models.CharField(maxlength=256) f_table_catalog = models.CharField(maxlength=256)
f_table_schema = models.CharField(maxlength=256) f_table_schema = models.CharField(maxlength=256)
f_table_name = models.CharField(maxlength=256) f_table_name = models.CharField(maxlength=256, primary_key=True)
f_geometry_column = models.CharField(maxlength=256) f_geometry_column = models.CharField(maxlength=256)
coord_dimension = models.IntegerField() coord_dimension = models.IntegerField()
srid = models.IntegerField() srid = models.IntegerField()
@ -29,6 +32,9 @@ class GeometryColumns(models.Model):
class Meta: class Meta:
db_table = 'geometry_columns' db_table = 'geometry_columns'
def __str__(self):
return "%s.%s - %dD %s field (SRID: %d)" % (self.f_table_name, self.f_geometry_column, self.coord_dimension, self.type, self.srid)
# This is the global 'spatial_ref_sys' table from PostGIS. # This is the global 'spatial_ref_sys' table from PostGIS.
# See PostGIS Documentation at Ch. 4.2.1 # See PostGIS Documentation at Ch. 4.2.1
class SpatialRefSys(models.Model): class SpatialRefSys(models.Model):
@ -36,32 +42,117 @@ class SpatialRefSys(models.Model):
auth_name = models.CharField(maxlength=256) auth_name = models.CharField(maxlength=256)
auth_srid = models.IntegerField() auth_srid = models.IntegerField()
srtext = models.CharField(maxlength=2048) srtext = models.CharField(maxlength=2048)
proj4text = models.CharField(maxlength=2048) proj4 = models.CharField(maxlength=2048, db_column='proj4text')
class Meta: class Meta:
db_table = 'spatial_ref_sys' db_table = 'spatial_ref_sys'
def _cache_osr(self):
"Caches a GDAL OSR object for this Spatial Reference."
if HAS_OSR:
if not hasattr(self, '_srs'):
# Trying to get from WKT first
try:
self._srs = SpatialReference(self.srtext, 'wkt')
return
except:
pass
# Trying the proj4 text next
try:
self._srs = SpatialReference(self.proj4, 'proj4')
return
except:
pass
raise Exception, 'Could not get a OSR Spatial Reference.'
else:
raise Exception, 'GDAL is not installed!'
@property
def srs(self):
self._cache_osr()
return self._srs.clone()
@property
def ellipsoid(self):
"""Returns a tuple of the ellipsoid parameters:
(semimajor axis, semiminor axis, and inverse flattening)."""
if HAS_OSR:
# Setting values initially to False
self._cache_osr()
major = self._srs.semi_major
minor = self._srs.semi_minor
invflat = self._srs.inverse_flattening
return (major, minor, invflat)
else:
m = spheroid_regex.match(self.srtext)
if m: return (float(m.group('major')), float(m.group('flattening')))
else: return None
@property
def name(self):
"Returns the projection name."
self._cache_osr()
return self._srs.name
@property @property
def spheroid(self): def spheroid(self):
"Pulls out the spheroid from the srtext." "Returns the spheroid for this spatial reference."
m = spheroid_regex.match(self.srtext) self._cache_osr()
if m: return self._srs['spheroid']
return (m.group('name'), float(m.group('major')), float(m.group('flattening')))
else:
return None
@property @property
def projected_units(self): def datum(self):
"If the spatial reference system is projected, get the units or return None." "Returns the datum for this spatial reference."
m = unit_regex.match(self.srtext) self._cache_osr()
if m: return self._srs['datum']
if m.group('authority'):
authority = (m.group('authority'), int(m.group('code'))) @property
else: def projected(self):
authority = None "Is this Spatial Reference projected?"
return (m.group('units'), float(m.group('conversion')), authority) self._cache_osr()
else: return self._srs.projected
return None
@property
def local(self):
"Is this Spatial Reference local?"
self._cache_osr()
return self._srs.local
@property
def geographic(self):
"Is this Spatial Reference geographic?"
self._cache_osr()
return self._srs.geographic
@property
def linear_name(self):
"Returns the linear units name."
self._cache_osr()
return self._srs.linear_name
@property
def linear_units(self):
"Returns the linear units."
self._cache_osr()
return self._srs.linear_units
@property
def angular_units(self):
"Returns the angular units."
self._cache_osr()
return self._srs.angular_units
@property
def angular_name(self):
"Returns the name of the angular units."
self._cache_osr()
return self._srs.angular_name
def __str__(self): def __str__(self):
return "%d - %s " % (self.srid, self.auth_name) "Returns the string representation. If GDAL is installed, it will be 'pretty' OGC WKT."
if HAS_OSR:
self._cache_osr()
if hasattr(self, '_srs'): return str(self._srs)
return "%d:%s " % (self.srid, self.auth_name)

View File

@ -0,0 +1,14 @@
from unittest import TestSuite, makeSuite, TextTestRunner
import test_geos, test_gdal_ds, test_gdal_srs, test_gdal_geom, test_spatialrefsys
def suite():
s = TestSuite()
s.addTest(test_geos.suite())
s.addTest(test_gdal_ds.suite())
s.addTest(test_gdal_srs.suite())
s.addTest(test_gdal_geom.suite())
s.addTest(test_spatialrefsys.suite())
return s
def run(verbosity=2):
TextTestRunner(verbosity=verbosity).run(suite())

View File

@ -0,0 +1,125 @@
import os, os.path, unittest
from django.contrib.gis.gdal import DataSource, OGRException
from django.contrib.gis.gdal.Field import OFTReal, OFTInteger, OFTString
# Path for SHP files
shp_path = os.path.dirname(__file__)
def get_shp(name):
return shp_path + os.sep + name + os.sep + name + '.shp'
# Test SHP data source object
class TestSHP:
def __init__(self, shp, **kwargs):
self.ds = get_shp(shp)
for key, value in kwargs.items():
setattr(self, key, value)
# List of acceptable data sources.
ds_list = (TestSHP('test_point', nfeat=5, nfld=3, geom='POINT', gtype=1, fields={'dbl' : OFTReal, 'int' : OFTReal, 'str' : OFTString,},
srs_wkt='GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'),
TestSHP('test_poly', nfeat=3, nfld=3, geom='POLYGON', gtype=3, fields={'float' : OFTReal, 'int' : OFTReal, 'str' : OFTString,},
srs_wkt='GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'),
)
bad_ds = (TestSHP('foo'),
)
class DataSourceTest(unittest.TestCase):
def test01_valid_shp(self):
"Testing valid SHP Data Source files."
for source in ds_list:
# Loading up the data source
ds = DataSource(source.ds)
# Making sure the layer count is what's expected (only 1 layer in a SHP file)
self.assertEqual(1, len(ds))
# Making sure GetName works
self.assertEqual(source.ds, ds.name)
# Making sure the driver name matches up
self.assertEqual('ESRI Shapefile', ds.driver)
# Making sure indexing works
try:
ds[len(ds)]
except IndexError:
pass
else:
self.fail('Expected an IndexError!')
def test02_invalid_shp(self):
"Testing invalid SHP files for the Data Source."
for source in bad_ds:
self.assertRaises(OGRException, DataSource, source.ds)
def test03_layers(self):
"Testing Data Source Layers."
for source in ds_list:
ds = DataSource(source.ds)
# Incrementing through each layer, this tests __iter__
for layer in ds:
# Making sure we get the number of features we expect
self.assertEqual(len(layer), source.nfeat)
# Making sure we get the number of fields we expect
self.assertEqual(layer.num_fields, source.nfld)
self.assertEqual(len(layer.fields), source.nfld)
# Now checking the field names.
flds = layer.fields
for f in flds: self.assertEqual(True, f in source.fields)
def test04_features(self):
"Testing Data Source Features."
for source in ds_list:
ds = DataSource(source.ds)
# Incrementing through each layer
for layer in ds:
# Incrementing through each feature in the layer
for feat in layer:
# Making sure the number of fields is what's expected.
self.assertEqual(source.nfld, len(feat))
self.assertEqual(source.gtype, feat.geom_type)
# Making sure the fields match to an appropriate OFT type.
for k, v in source.fields.items():
fld = feat[k] # Indexing with string value
# Asserting the string representation (which asserts the type)
self.assertEqual('%s (%s)' % (k, v.__name__), str(fld))
# Testing __iter__ on the Feature
for fld in feat: self.assertEqual(fld.name in source.fields.keys(), True)
def test05_geometries(self):
"Testing Geometries from Data Source Features."
for source in ds_list:
ds = DataSource(source.ds)
# Incrementing through each layer and feature.
for layer in ds:
for feat in layer:
g = feat.geom
# Making sure we get the right Geometry name & type
self.assertEqual(source.geom, g.geom_name)
self.assertEqual(source.gtype, g.geom_type)
# Making sure the SpatialReference is as expected.
self.assertEqual(source.srs_wkt, g.srs.wkt)
def suite():
s = unittest.TestSuite()
s.addTest(unittest.makeSuite(DataSourceTest))
return s
def run(verbosity=2):
unittest.TextTestRunner(verbosity=verbosity).run(suite())

View File

@ -0,0 +1,79 @@
import unittest
from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, OGRException
from geometries import *
class OGRGeomTest(unittest.TestCase):
"This tests the OGR Geometry."
def test00_geomtype(self):
"Testing OGRGeomType object."
# OGRGeomType should initialize on all these inputs.
try:
g = OGRGeomType(0)
g = OGRGeomType(1)
g = OGRGeomType(7)
g = OGRGeomType('point')
g = OGRGeomType('GeometrycollectioN')
except:
self.fail('Could not create an OGRGeomType object!')
# Should throw TypeError on this input
self.assertRaises(TypeError, OGRGeomType.__init__, 23)
self.assertRaises(TypeError, OGRGeomType.__init__, 'fooD')
self.assertRaises(TypeError, OGRGeomType.__init__, 9)
# Equivalence can take strings, ints, and other OGRGeomTypes
self.assertEqual(True, OGRGeomType(1) == OGRGeomType(1))
self.assertEqual(True, OGRGeomType(7) == 'GeometryCollection')
self.assertEqual(True, OGRGeomType('point') == 'POINT')
self.assertEqual(False, OGRGeomType('point') == 2)
self.assertEqual(True, OGRGeomType(6) == 'MULtiPolyGON')
def test01_wkt(self):
"Testing WKT output."
for g in wkt_out:
geom = OGRGeometry(g.wkt)
def test02_points(self):
"Testing Point objects."
prev = OGRGeometry('POINT(0 0)')
for p in points:
if not hasattr(p, 'z'): # No 3D
pnt = OGRGeometry(p.wkt)
self.assertEqual(pnt.geom_type, 1)
self.assertEqual(p.x, pnt.x)
self.assertEqual(p.y, pnt.y)
self.assertEqual((p.x, p.y), pnt.tuple)
def test03_polygons(self):
"Testing Polygon objects."
for p in polygons:
poly = OGRGeometry(p.wkt)
first = True
for r in poly:
if first and p.ext_ring_cs:
first = False
# Testing the equivilance of the exerior rings
# since the first iteration will be the exterior ring.
self.assertEqual(len(p.ext_ring_cs), r.point_count)
self.assertEqual(p.ext_ring_cs, r.tuple)
def test04_multipoints(self):
"Testing MultiPoint objects."
for mp in multipoints:
mgeom1 = OGRGeometry(mp.wkt) # First one from WKT
mgeom2 = OGRGeometry('MULTIPOINT') # Creating empty multipoint
for g in mgeom1:
mgeom2.add(g) # adding each point from the multipoint
self.assertEqual(mgeom1, mgeom2) # they should equal
def suite():
s = unittest.TestSuite()
s.addTest(unittest.makeSuite(OGRGeomTest))
return s
def run(verbosity=2):
unittest.TextTestRunner(verbosity=verbosity).run(suite())

View File

@ -0,0 +1,144 @@
import unittest
from django.contrib.gis.gdal import SpatialReference, CoordTransform, OGRException, SRSException
class TestSRS:
def __init__(self, wkt, **kwargs):
self.wkt = wkt
for key, value in kwargs.items():
setattr(self, key, value)
# Some Spatial Reference examples
srlist = (TestSRS('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]',
proj='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ',
epsg=4326, projected=False, geographic=True, local=False,
lin_name='unknown', ang_name='degree', lin_units=1.0, ang_units=0.0174532925199,
auth={'GEOGCS' : ('EPSG', '4326'), 'spheroid' : ('EPSG', '7030')},
attr=(('DATUM', 'WGS_1984'), (('SPHEROID', 1), '6378137'),('primem|authority', 'EPSG'),),
),
TestSRS('PROJCS["NAD83 / Texas South Central",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",30.28333333333333],PARAMETER["standard_parallel_2",28.38333333333333],PARAMETER["latitude_of_origin",27.83333333333333],PARAMETER["central_meridian",-99],PARAMETER["false_easting",600000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32140"]]',
proj='+proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ',
epsg=32140, projected=True, geographic=False, local=False,
lin_name='metre', ang_name='degree', lin_units=1.0, ang_units=0.0174532925199,
auth={'PROJCS' : ('EPSG', '32140'), 'spheroid' : ('EPSG', '7019'), 'unit' : ('EPSG', '9001'),},
attr=(('DATUM', 'North_American_Datum_1983'),(('SPHEROID', 2), '298.257222101'),('PROJECTION','Lambert_Conformal_Conic_2SP'),),
),
TestSRS('PROJCS["NAD_1983_StatePlane_Texas_South_Central_FIPS_4204_Feet",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",1968500.0],PARAMETER["False_Northing",13123333.33333333],PARAMETER["Central_Meridian",-99.0],PARAMETER["Standard_Parallel_1",28.38333333333333],PARAMETER["Standard_Parallel_2",30.28333333333334],PARAMETER["Latitude_Of_Origin",27.83333333333333],UNIT["Foot_US",0.3048006096012192]]',
proj='+proj=lcc +lat_1=28.38333333333333 +lat_2=30.28333333333334 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=3999999.999999999 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs ',
epsg=None, projected=True, geographic=False, local=False,
lin_name='Foot_US', ang_name='Degree', lin_units=0.3048006096012192, ang_units=0.0174532925199,
auth={'PROJCS' : (None, None),},
attr=(('PROJCS|GeOgCs|spheroid', 'GRS_1980'),(('projcs', 9), 'UNIT'), (('projcs', 11), None),),
),
# This is really ESRI format, not WKT -- but the import should work the same
TestSRS('LOCAL_CS["Non-Earth (Meter)",LOCAL_DATUM["Local Datum",0],UNIT["Meter",1.0],AXIS["X",EAST],AXIS["Y",NORTH]]',
esri=True, proj=None, epsg=None, projected=False, geographic=False, local=True,
lin_name='Meter', ang_name='degree', lin_units=1.0, ang_units=0.0174532925199,
attr=(('LOCAL_DATUM', 'Local Datum'), ('unit', 'Meter')),
),
)
# Well-Known Names
well_known = (TestSRS('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]', wk='WGS84'),
TestSRS('GEOGCS["WGS 72",DATUM["WGS_1972",SPHEROID["WGS 72",6378135,298.26,AUTHORITY["EPSG","7043"]],AUTHORITY["EPSG","6322"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4322"]]', wk='WGS72'),
TestSRS('GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]]', wk='NAD27'),
TestSRS('GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]]', wk='NAD83'),
TestSRS('PROJCS["NZGD49 / Karamea Circuit",GEOGCS["NZGD49",DATUM["New_Zealand_Geodetic_Datum_1949",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993],AUTHORITY["EPSG","6272"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4272"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",-41.28991152777778],PARAMETER["central_meridian",172.1090281944444],PARAMETER["scale_factor",1],PARAMETER["false_easting",300000],PARAMETER["false_northing",700000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","27216"]]', wk='EPSG:27216'),
)
bad_srlist = ('Foobar', 'OOJCS["NAD83 / Texas South Central",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",30.28333333333333],PARAMETER["standard_parallel_2",28.38333333333333],PARAMETER["latitude_of_origin",27.83333333333333],PARAMETER["central_meridian",-99],PARAMETER["false_easting",600000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32140"]]',)
class SpatialRefTest(unittest.TestCase):
def test01_wkt(self):
"Testing initialization on valid OGC WKT."
for s in srlist:
srs = SpatialReference(s.wkt)
def test02_bad_wkt(self):
"Testing initialization on invalid WKT."
for bad in bad_srlist:
try:
srs = SpatialReference(bad)
srs.validate()
except (SRSException, OGRException):
pass
else:
self.fail('Should not have initialized on bad WKT "%s"!')
def test03_get_wkt(self):
"Testing getting the WKT."
for s in srlist:
srs = SpatialReference(s.wkt)
self.assertEqual(s.wkt, srs.wkt)
def test04_proj(self):
"Test PROJ.4 import and export."
for s in srlist:
if s.proj:
srs1 = SpatialReference(s.wkt)
srs2 = SpatialReference(s.proj, 'proj')
self.assertEqual(srs1.proj, srs2.proj)
def test05_epsg(self):
"Test EPSG import."
for s in srlist:
if s.epsg:
srs1 = SpatialReference(s.wkt)
srs2 = SpatialReference(s.epsg, 'epsg')
self.assertEqual(srs1.wkt, srs2.wkt)
def test07_boolean_props(self):
"Testing the boolean properties."
for s in srlist:
srs = SpatialReference(s.wkt)
self.assertEqual(s.projected, srs.projected)
self.assertEqual(s.geographic, srs.geographic)
def test08_angular_linear(self):
"Testing the linear and angular units routines."
for s in srlist:
srs = SpatialReference(s.wkt)
self.assertEqual(s.ang_name, srs.angular_name)
self.assertEqual(s.lin_name, srs.linear_name)
self.assertAlmostEqual(s.ang_units, srs.angular_units, 9)
self.assertAlmostEqual(s.lin_units, srs.linear_units, 9)
def test09_authority(self):
"Testing the authority name & code routines."
for s in srlist:
if hasattr(s, 'auth'):
srs = SpatialReference(s.wkt)
for target, tup in s.auth.items():
self.assertEqual(tup[0], srs.auth_name(target))
self.assertEqual(tup[1], srs.auth_code(target))
def test10_attributes(self):
"Testing the attribute retrieval routines."
for s in srlist:
srs = SpatialReference(s.wkt)
for tup in s.attr:
att = tup[0] # Attribute to test
exp = tup[1] # Expected result
self.assertEqual(exp, srs[att])
def test11_wellknown(self):
"Testing Well Known Names of Spatial References."
for s in well_known:
srs = SpatialReference(s.wk)
self.assertEqual(s.wkt, srs.wkt)
def test12_coordtransform(self):
"Testing initialization of a CoordTransform."
target = SpatialReference('WGS84')
for s in srlist:
if s.proj:
ct = CoordTransform(SpatialReference(s.wkt), target)
def suite():
s = unittest.TestSuite()
s.addTest(unittest.makeSuite(SpatialRefTest))
return s
def run(verbosity=2):
unittest.TextTestRunner(verbosity=verbosity).run(suite())

Binary file not shown.

View File

@ -0,0 +1 @@
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1 @@
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1,74 @@
import unittest
from django.contrib.gis.models import SpatialRefSys
test_srs = ({'srid' : 4326,
'auth_name' : 'EPSG',
'auth_srid' : 4326,
'srtext' : 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]',
'proj4' : '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ',
'spheroid' : 'WGS 84', 'name' : 'WGS 84',
'geographic' : True, 'projected' : False,
'ellipsoid' : (6378137.0, 6356752.3, 298.257223563), # From proj's "cs2cs -le" and Wikipedia (semi-minor only)
'eprec' : (1, 1, 9),
},
{'srid' : 32140,
'auth_name' : 'EPSG',
'auth_srid' : 32140,
'srtext' : 'PROJCS["NAD83 / Texas South Central",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",30.28333333333333],PARAMETER["standard_parallel_2",28.38333333333333],PARAMETER["latitude_of_origin",27.83333333333333],PARAMETER["central_meridian",-99],PARAMETER["false_easting",600000],PARAMETER["false_northing",4000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32140"]]',
'proj4' : '+proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ',
'spheroid' : 'GRS 1980', 'name' : 'NAD83 / Texas South Central',
'geographic' : False, 'projected' : True,
'ellipsoid' : (6378137.0, 6356752.31414, 298.257222101), # From proj's "cs2cs -le" and Wikipedia (semi-minor only)
'eprec' : (1, 5, 10),
},
)
class SpatialRefSysTest(unittest.TestCase):
def test01_retrieve(self):
"Testing retrieval of SpatialRefSys model objects."
for sd in test_srs:
srs = SpatialRefSys.objects.get(srid=sd['srid'])
self.assertEqual(srs.srid, sd['srid'])
self.assertEqual(srs.auth_name, sd['auth_name'])
self.assertEqual(srs.auth_srid, sd['auth_srid'])
self.assertEqual(srs.srtext, sd['srtext'])
self.assertEqual(srs.proj4, sd['proj4'])
def test02_osr(self):
"Testing getting OSR objects from SpatialRefSys model objects."
for sd in test_srs:
sr = SpatialRefSys.objects.get(srid=sd['srid'])
self.assertEqual(sd['spheroid'], sr.spheroid)
self.assertEqual(sd['geographic'], sr.geographic)
self.assertEqual(sd['projected'], sr.projected)
self.assertEqual(sd['name'], sr.name)
# Testing the SpatialReference object directly.
srs = sr.srs
self.assertEqual(sd['proj4'], srs.proj)
self.assertEqual(sd['srtext'], srs.wkt)
def test03_ellipsoid(self):
"Testing the ellipsoid property."
for sd in test_srs:
# Getting the ellipsoid and precision parameters.
ellps1 = sd['ellipsoid']
prec = sd['eprec']
# Getting our spatial reference and its ellipsoid
srs = SpatialRefSys.objects.get(srid=sd['srid'])
ellps2 = srs.ellipsoid
for i in range(3):
param1 = ellps1[i]
param2 = ellps2[i]
self.assertAlmostEqual(ellps1[i], ellps2[i], prec[i])
def suite():
s = unittest.TestSuite()
s.addTest(unittest.makeSuite(SpatialRefSysTest))
return s
def run(verbosity=2):
unittest.TextTestRunner(verbosity=verbosity).run(suite())