before send to remote

This commit is contained in:
2022-08-26 16:41:18 +06:00
commit 3814beb3e0
5408 changed files with 652023 additions and 0 deletions

View File

@@ -0,0 +1,28 @@
Copyright (c) 2007-2009, 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,58 @@
"""
This module houses ctypes interfaces for GDAL objects. The following GDAL
objects are supported:
CoordTransform: Used for coordinate transformations from one spatial
reference system to another.
Driver: Wraps an OGR data source driver.
DataSource: Wrapper for the OGR data source object, supports
OGR-supported data sources.
Envelope: A ctypes structure for bounding boxes (GDAL library
not required).
OGRGeometry: Object for accessing OGR Geometry functionality.
OGRGeomType: A class for representing the different OGR Geometry
types (GDAL library not required).
SpatialReference: Represents OSR Spatial Reference objects.
The GDAL library will be imported from the system path using the default
library name for the current OS. The default library path may be overridden
by setting `GDAL_LIBRARY_PATH` in your settings with the path to the GDAL C
library on your system.
"""
from django.contrib.gis.gdal.datasource import DataSource
from django.contrib.gis.gdal.driver import Driver
from django.contrib.gis.gdal.envelope import Envelope
from django.contrib.gis.gdal.error import GDALException, SRSException, check_err
from django.contrib.gis.gdal.geometries import OGRGeometry
from django.contrib.gis.gdal.geomtype import OGRGeomType
from django.contrib.gis.gdal.libgdal import (
GDAL_VERSION,
gdal_full_version,
gdal_version,
)
from django.contrib.gis.gdal.raster.source import GDALRaster
from django.contrib.gis.gdal.srs import AxisOrder, CoordTransform, SpatialReference
__all__ = (
"AxisOrder",
"Driver",
"DataSource",
"CoordTransform",
"Envelope",
"GDALException",
"GDALRaster",
"GDAL_VERSION",
"OGRGeometry",
"OGRGeomType",
"SpatialReference",
"SRSException",
"check_err",
"gdal_version",
"gdal_full_version",
)

View File

@@ -0,0 +1,6 @@
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.ptr import CPointerBase
class GDALBase(CPointerBase):
null_ptr_exception_class = GDALException

View File

@@ -0,0 +1,126 @@
"""
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).
When instantiating a DataSource object, use the filename of a
GDAL-supported data source. For example, a SHP file or a
TIGER/Line file from the government.
The ds_driver keyword is used internally when a ctypes pointer
is passed in directly.
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']
# We can also increment through all of the fields
# attached to this feature.
for field in feature:
# Get the name of the field (e.g. 'description')
nm = field.name
# Get the type (integer) of the field, e.g. 0 => OFTInteger
t = field.type
# Returns the value the field; OFTIntegers return ints,
# OFTReal returns floats, all else returns string.
val = field.value
"""
from ctypes import byref
from pathlib import Path
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.driver import Driver
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.layer import Layer
from django.contrib.gis.gdal.prototypes import ds as capi
from django.utils.encoding import force_bytes, force_str
# For more information, see the OGR C API documentation:
# https://gdal.org/api/vector_c_api.html
#
# The OGR_DS_* routines are relevant here.
class DataSource(GDALBase):
"Wraps an OGR Data Source object."
destructor = capi.destroy_ds
def __init__(self, ds_input, ds_driver=False, write=False, encoding="utf-8"):
# The write flag.
if write:
self._write = 1
else:
self._write = 0
# See also https://gdal.org/development/rfc/rfc23_ogr_unicode.html
self.encoding = encoding
Driver.ensure_registered()
if isinstance(ds_input, (str, Path)):
# The data source driver is a void pointer.
ds_driver = Driver.ptr_type()
try:
# OGROpen will auto-detect the data source type.
ds = capi.open_ds(force_bytes(ds_input), self._write, byref(ds_driver))
except GDALException:
# Making the error message more clear rather than something
# like "Invalid pointer returned from OGROpen".
raise GDALException('Could not open the datasource at "%s"' % ds_input)
elif isinstance(ds_input, self.ptr_type) and isinstance(
ds_driver, Driver.ptr_type
):
ds = ds_input
else:
raise GDALException("Invalid data source input type: %s" % type(ds_input))
if ds:
self.ptr = ds
self.driver = Driver(ds_driver)
else:
# Raise an exception if the returned pointer is NULL
raise GDALException('Invalid data source file "%s"' % ds_input)
def __getitem__(self, index):
"Allows use of the index [] operator to get a layer at the index."
if isinstance(index, str):
try:
layer = capi.get_layer_by_name(self.ptr, force_bytes(index))
except GDALException:
raise IndexError("Invalid OGR layer name given: %s." % index)
elif isinstance(index, int):
if 0 <= index < self.layer_count:
layer = capi.get_layer(self._ptr, index)
else:
raise IndexError(
"Index out of range when accessing layers in a datasource: %s."
% index
)
else:
raise TypeError("Invalid index type: %s" % type(index))
return Layer(layer, self)
def __len__(self):
"Return the number of layers within the data source."
return self.layer_count
def __str__(self):
"Return OGR GetName and Driver for the Data Source."
return "%s (%s)" % (self.name, self.driver)
@property
def layer_count(self):
"Return the number of layers in the data source."
return capi.get_layer_count(self._ptr)
@property
def name(self):
"Return the name of the data source."
name = capi.get_ds_name(self._ptr)
return force_str(name, self.encoding, strings_only=True)

View File

@@ -0,0 +1,103 @@
from ctypes import c_void_p
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.prototypes import ds as vcapi
from django.contrib.gis.gdal.prototypes import raster as rcapi
from django.utils.encoding import force_bytes, force_str
class Driver(GDALBase):
"""
Wrap a GDAL/OGR Data Source Driver.
For more information, see the C API documentation:
https://gdal.org/api/vector_c_api.html
https://gdal.org/api/raster_c_api.html
"""
# Case-insensitive aliases for some GDAL/OGR Drivers.
# For a complete list of original driver names see
# https://gdal.org/drivers/vector/
# https://gdal.org/drivers/raster/
_alias = {
# vector
"esri": "ESRI Shapefile",
"shp": "ESRI Shapefile",
"shape": "ESRI Shapefile",
"tiger": "TIGER",
"tiger/line": "TIGER",
# raster
"tiff": "GTiff",
"tif": "GTiff",
"jpeg": "JPEG",
"jpg": "JPEG",
}
def __init__(self, dr_input):
"""
Initialize an GDAL/OGR driver on either a string or integer input.
"""
if isinstance(dr_input, str):
# If a string name of the driver was passed in
self.ensure_registered()
# Checking the alias dictionary (case-insensitive) to see if an
# alias exists for the given driver.
if dr_input.lower() in self._alias:
name = self._alias[dr_input.lower()]
else:
name = dr_input
# Attempting to get the GDAL/OGR driver by the string name.
for iface in (vcapi, rcapi):
driver = c_void_p(iface.get_driver_by_name(force_bytes(name)))
if driver:
break
elif isinstance(dr_input, int):
self.ensure_registered()
for iface in (vcapi, rcapi):
driver = iface.get_driver(dr_input)
if driver:
break
elif isinstance(dr_input, c_void_p):
driver = dr_input
else:
raise GDALException(
"Unrecognized input type for GDAL/OGR Driver: %s" % type(dr_input)
)
# Making sure we get a valid pointer to the OGR Driver
if not driver:
raise GDALException(
"Could not initialize GDAL/OGR Driver on input: %s" % dr_input
)
self.ptr = driver
def __str__(self):
return self.name
@classmethod
def ensure_registered(cls):
"""
Attempt to register all the data source drivers.
"""
# Only register all if the driver counts are 0 (or else all drivers
# will be registered over and over again)
if not vcapi.get_driver_count():
vcapi.register_all()
if not rcapi.get_driver_count():
rcapi.register_all()
@classmethod
def driver_count(cls):
"""
Return the number of GDAL/OGR data source drivers registered.
"""
return vcapi.get_driver_count() + rcapi.get_driver_count()
@property
def name(self):
"""
Return description/name string for this driver.
"""
return force_str(rcapi.get_driver_description(self.ptr))

View File

@@ -0,0 +1,203 @@
"""
The GDAL/OGR library uses an Envelope structure to hold the bounding
box information for a geometry. The envelope (bounding box) contains
two pairs of coordinates, one for the lower left coordinate and one
for the upper right coordinate:
+----------o Upper right; (max_x, max_y)
| |
| |
| |
Lower left (min_x, min_y) o----------+
"""
from ctypes import Structure, c_double
from django.contrib.gis.gdal.error import GDALException
# The OGR definition of an Envelope is a C structure containing four doubles.
# See the 'ogr_core.h' source file for more information:
# https://gdal.org/doxygen/ogr__core_8h_source.html
class OGREnvelope(Structure):
"Represent the OGREnvelope C Structure."
_fields_ = [
("MinX", c_double),
("MaxX", c_double),
("MinY", c_double),
("MaxY", c_double),
]
class Envelope:
"""
The Envelope object is a C structure that contains the minimum and
maximum X, Y coordinates for a rectangle bounding box. The naming
of the variables is compatible with the OGR Envelope structure.
"""
def __init__(self, *args):
"""
The initialization function may take an OGREnvelope structure, 4-element
tuple or list, or 4 individual arguments.
"""
if len(args) == 1:
if isinstance(args[0], OGREnvelope):
# OGREnvelope (a ctypes Structure) was passed in.
self._envelope = args[0]
elif isinstance(args[0], (tuple, list)):
# A tuple was passed in.
if len(args[0]) != 4:
raise GDALException(
"Incorrect number of tuple elements (%d)." % len(args[0])
)
else:
self._from_sequence(args[0])
else:
raise TypeError("Incorrect type of argument: %s" % type(args[0]))
elif len(args) == 4:
# Individual parameters passed in.
# Thanks to ww for the help
self._from_sequence([float(a) for a in args])
else:
raise GDALException("Incorrect number (%d) of arguments." % len(args))
# Checking the x,y coordinates
if self.min_x > self.max_x:
raise GDALException("Envelope minimum X > maximum X.")
if self.min_y > self.max_y:
raise GDALException("Envelope minimum Y > maximum Y.")
def __eq__(self, other):
"""
Return True if the envelopes are equivalent; can compare against
other Envelopes and 4-tuples.
"""
if isinstance(other, Envelope):
return (
(self.min_x == other.min_x)
and (self.min_y == other.min_y)
and (self.max_x == other.max_x)
and (self.max_y == other.max_y)
)
elif isinstance(other, tuple) and len(other) == 4:
return (
(self.min_x == other[0])
and (self.min_y == other[1])
and (self.max_x == other[2])
and (self.max_y == other[3])
)
else:
raise GDALException("Equivalence testing only works with other Envelopes.")
def __str__(self):
"Return a string representation of the tuple."
return str(self.tuple)
def _from_sequence(self, seq):
"Initialize the C OGR Envelope structure from the given sequence."
self._envelope = OGREnvelope()
self._envelope.MinX = seq[0]
self._envelope.MinY = seq[1]
self._envelope.MaxX = seq[2]
self._envelope.MaxY = seq[3]
def expand_to_include(self, *args):
"""
Modify the envelope to expand to include the boundaries of
the passed-in 2-tuple (a point), 4-tuple (an extent) or
envelope.
"""
# We provide a number of different signatures for this method,
# and the logic here is all about converting them into a
# 4-tuple single parameter which does the actual work of
# expanding the envelope.
if len(args) == 1:
if isinstance(args[0], Envelope):
return self.expand_to_include(args[0].tuple)
elif hasattr(args[0], "x") and hasattr(args[0], "y"):
return self.expand_to_include(
args[0].x, args[0].y, args[0].x, args[0].y
)
elif isinstance(args[0], (tuple, list)):
# A tuple was passed in.
if len(args[0]) == 2:
return self.expand_to_include(
(args[0][0], args[0][1], args[0][0], args[0][1])
)
elif len(args[0]) == 4:
(minx, miny, maxx, maxy) = args[0]
if minx < self._envelope.MinX:
self._envelope.MinX = minx
if miny < self._envelope.MinY:
self._envelope.MinY = miny
if maxx > self._envelope.MaxX:
self._envelope.MaxX = maxx
if maxy > self._envelope.MaxY:
self._envelope.MaxY = maxy
else:
raise GDALException(
"Incorrect number of tuple elements (%d)." % len(args[0])
)
else:
raise TypeError("Incorrect type of argument: %s" % type(args[0]))
elif len(args) == 2:
# An x and an y parameter were passed in
return self.expand_to_include((args[0], args[1], args[0], args[1]))
elif len(args) == 4:
# Individual parameters passed in.
return self.expand_to_include(args)
else:
raise GDALException("Incorrect number (%d) of arguments." % len(args[0]))
@property
def min_x(self):
"Return the value of the minimum X coordinate."
return self._envelope.MinX
@property
def min_y(self):
"Return the value of the minimum Y coordinate."
return self._envelope.MinY
@property
def max_x(self):
"Return the value of the maximum X coordinate."
return self._envelope.MaxX
@property
def max_y(self):
"Return the value of the maximum Y coordinate."
return self._envelope.MaxY
@property
def ur(self):
"Return the upper-right coordinate."
return (self.max_x, self.max_y)
@property
def ll(self):
"Return the lower-left coordinate."
return (self.min_x, self.min_y)
@property
def tuple(self):
"Return a tuple representing the envelope."
return (self.min_x, self.min_y, self.max_x, self.max_y)
@property
def wkt(self):
"Return WKT representing a Polygon for this envelope."
# TODO: Fix significant figures.
return "POLYGON((%s %s,%s %s,%s %s,%s %s,%s %s))" % (
self.min_x,
self.min_y,
self.min_x,
self.max_y,
self.max_x,
self.max_y,
self.max_x,
self.min_y,
self.min_x,
self.min_y,
)

View File

@@ -0,0 +1,61 @@
"""
This module houses the GDAL & SRS Exception objects, and the
check_err() routine which checks the status code returned by
GDAL/OGR methods.
"""
# #### GDAL & SRS Exceptions ####
class GDALException(Exception):
pass
class SRSException(Exception):
pass
# #### GDAL/OGR error checking codes and routine ####
# OGR Error Codes
OGRERR_DICT = {
1: (GDALException, "Not enough data."),
2: (GDALException, "Not enough memory."),
3: (GDALException, "Unsupported geometry type."),
4: (GDALException, "Unsupported operation."),
5: (GDALException, "Corrupt data."),
6: (GDALException, "OGR failure."),
7: (SRSException, "Unsupported SRS."),
8: (GDALException, "Invalid handle."),
}
# CPL Error Codes
# https://gdal.org/api/cpl.html#cpl-error-h
CPLERR_DICT = {
1: (GDALException, "AppDefined"),
2: (GDALException, "OutOfMemory"),
3: (GDALException, "FileIO"),
4: (GDALException, "OpenFailed"),
5: (GDALException, "IllegalArg"),
6: (GDALException, "NotSupported"),
7: (GDALException, "AssertionFailed"),
8: (GDALException, "NoWriteAccess"),
9: (GDALException, "UserInterrupt"),
10: (GDALException, "ObjectNull"),
}
ERR_NONE = 0
def check_err(code, cpl=False):
"""
Check the given CPL/OGRERR and raise an exception where appropriate.
"""
err_dict = CPLERR_DICT if cpl else OGRERR_DICT
if code == ERR_NONE:
return
elif code in err_dict:
e, msg = err_dict[code]
raise e(msg)
else:
raise GDALException('Unknown error code: "%s"' % code)

View File

@@ -0,0 +1,120 @@
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.field import Field
from django.contrib.gis.gdal.geometries import OGRGeometry, OGRGeomType
from django.contrib.gis.gdal.prototypes import ds as capi
from django.contrib.gis.gdal.prototypes import geom as geom_api
from django.utils.encoding import force_bytes, force_str
# For more information, see the OGR C API source code:
# https://gdal.org/api/vector_c_api.html
#
# The OGR_F_* routines are relevant here.
class Feature(GDALBase):
"""
This class that wraps an OGR Feature, needs to be instantiated
from a Layer object.
"""
destructor = capi.destroy_feature
def __init__(self, feat, layer):
"""
Initialize Feature from a pointer and its Layer object.
"""
if not feat:
raise GDALException("Cannot create OGR Feature, invalid pointer given.")
self.ptr = feat
self._layer = layer
def __getitem__(self, index):
"""
Get the Field object at the specified index, which may be either
an integer or the Field's string label. Note that the Field object
is not the field's _value_ -- use the `get` method instead to
retrieve the value (e.g. an integer) instead of a Field instance.
"""
if isinstance(index, str):
i = self.index(index)
elif 0 <= index < self.num_fields:
i = index
else:
raise IndexError(
"Index out of range when accessing field in a feature: %s." % index
)
return Field(self, i)
def __len__(self):
"Return 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):
"Do equivalence testing on the features."
return bool(capi.feature_equal(self.ptr, other._ptr))
# #### Feature Properties ####
@property
def encoding(self):
return self._layer._ds.encoding
@property
def fid(self):
"Return the feature identifier."
return capi.get_fid(self.ptr)
@property
def layer_name(self):
"Return the name of the layer for the feature."
name = capi.get_feat_name(self._layer._ldefn)
return force_str(name, self.encoding, strings_only=True)
@property
def num_fields(self):
"Return the number of fields in the Feature."
return capi.get_feat_field_count(self.ptr)
@property
def fields(self):
"Return a list of fields in the Feature."
return [
force_str(
capi.get_field_name(capi.get_field_defn(self._layer._ldefn, i)),
self.encoding,
strings_only=True,
)
for i in range(self.num_fields)
]
@property
def geom(self):
"Return the OGR Geometry for this Feature."
# Retrieving the geometry pointer for the feature.
geom_ptr = capi.get_feat_geom_ref(self.ptr)
return OGRGeometry(geom_api.clone_geom(geom_ptr))
@property
def geom_type(self):
"Return the OGR Geometry Type for this Feature."
return OGRGeomType(capi.get_fd_geom_type(self._layer._ldefn))
# #### Feature Methods ####
def get(self, field):
"""
Return the value of the field, instead of an instance of the Field
object. May take a string of the field name or a Field object as
parameters.
"""
field_name = getattr(field, "name", field)
return self[field_name].value
def index(self, field_name):
"Return the index of the given field name."
i = capi.get_field_index(self.ptr, force_bytes(field_name))
if i < 0:
raise IndexError("Invalid OFT field name given: %s." % field_name)
return i

View File

@@ -0,0 +1,253 @@
from ctypes import byref, c_int
from datetime import date, datetime, time
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.prototypes import ds as capi
from django.utils.encoding import force_str
# For more information, see the OGR C API source code:
# https://gdal.org/api/vector_c_api.html
#
# The OGR_Fld_* routines are relevant here.
class Field(GDALBase):
"""
Wrap an OGR Field. Needs to be instantiated from a Feature object.
"""
def __init__(self, feat, index):
"""
Initialize on the feature object and the integer index of
the field within the feature.
"""
# Setting the feature pointer and index.
self._feat = feat
self._index = index
# Getting the pointer for this field.
fld_ptr = capi.get_feat_field_defn(feat.ptr, index)
if not fld_ptr:
raise GDALException("Cannot create OGR Field, invalid pointer given.")
self.ptr = fld_ptr
# Setting the class depending upon the OGR Field Type (OFT)
self.__class__ = OGRFieldTypes[self.type]
def __str__(self):
"Return the string representation of the Field."
return str(self.value).strip()
# #### Field Methods ####
def as_double(self):
"Retrieve the Field's value as a double (float)."
return (
capi.get_field_as_double(self._feat.ptr, self._index)
if self.is_set
else None
)
def as_int(self, is_64=False):
"Retrieve the Field's value as an integer."
if is_64:
return (
capi.get_field_as_integer64(self._feat.ptr, self._index)
if self.is_set
else None
)
else:
return (
capi.get_field_as_integer(self._feat.ptr, self._index)
if self.is_set
else None
)
def as_string(self):
"Retrieve the Field's value as a string."
if not self.is_set:
return None
string = capi.get_field_as_string(self._feat.ptr, self._index)
return force_str(string, encoding=self._feat.encoding, strings_only=True)
def as_datetime(self):
"Retrieve the Field's value as a tuple of date & time components."
if not self.is_set:
return None
yy, mm, dd, hh, mn, ss, tz = [c_int() for i in range(7)]
status = capi.get_field_as_datetime(
self._feat.ptr,
self._index,
byref(yy),
byref(mm),
byref(dd),
byref(hh),
byref(mn),
byref(ss),
byref(tz),
)
if status:
return (yy, mm, dd, hh, mn, ss, tz)
else:
raise GDALException(
"Unable to retrieve date & time information from the field."
)
# #### Field Properties ####
@property
def is_set(self):
"Return True if the value of this field isn't null, False otherwise."
return capi.is_field_set(self._feat.ptr, self._index)
@property
def name(self):
"Return the name of this Field."
name = capi.get_field_name(self.ptr)
return force_str(name, encoding=self._feat.encoding, strings_only=True)
@property
def precision(self):
"Return the precision of this Field."
return capi.get_field_precision(self.ptr)
@property
def type(self):
"Return the OGR type of this Field."
return capi.get_field_type(self.ptr)
@property
def type_name(self):
"Return the OGR field type name for this Field."
return capi.get_field_type_name(self.type)
@property
def value(self):
"Return the value of this Field."
# Default is to get the field as a string.
return self.as_string()
@property
def width(self):
"Return the width of this Field."
return capi.get_field_width(self.ptr)
# ### The Field sub-classes for each OGR Field type. ###
class OFTInteger(Field):
_bit64 = False
@property
def value(self):
"Return an integer contained in this field."
return self.as_int(self._bit64)
@property
def type(self):
"""
GDAL uses OFTReals to represent OFTIntegers in created
shapefiles -- forcing the type here since the underlying field
type may actually be OFTReal.
"""
return 0
class OFTReal(Field):
@property
def value(self):
"Return a float contained in this field."
return self.as_double()
# String & Binary fields, just subclasses
class OFTString(Field):
pass
class OFTWideString(Field):
pass
class OFTBinary(Field):
pass
# OFTDate, OFTTime, OFTDateTime fields.
class OFTDate(Field):
@property
def value(self):
"Return a Python `date` object for the OFTDate field."
try:
yy, mm, dd, hh, mn, ss, tz = self.as_datetime()
return date(yy.value, mm.value, dd.value)
except (TypeError, ValueError, GDALException):
return None
class OFTDateTime(Field):
@property
def value(self):
"Return a Python `datetime` object for this OFTDateTime field."
# TODO: Adapt timezone information.
# See https://lists.osgeo.org/pipermail/gdal-dev/2006-February/007990.html
# The `tz` variable has values of: 0=unknown, 1=localtime (ambiguous),
# 100=GMT, 104=GMT+1, 80=GMT-5, etc.
try:
yy, mm, dd, hh, mn, ss, tz = self.as_datetime()
return datetime(yy.value, mm.value, dd.value, hh.value, mn.value, ss.value)
except (TypeError, ValueError, GDALException):
return None
class OFTTime(Field):
@property
def value(self):
"Return a Python `time` object for this OFTTime field."
try:
yy, mm, dd, hh, mn, ss, tz = self.as_datetime()
return time(hh.value, mn.value, ss.value)
except (ValueError, GDALException):
return None
class OFTInteger64(OFTInteger):
_bit64 = True
# List fields are also just subclasses
class OFTIntegerList(Field):
pass
class OFTRealList(Field):
pass
class OFTStringList(Field):
pass
class OFTWideStringList(Field):
pass
class OFTInteger64List(Field):
pass
# Class mapping dictionary for OFT Types and reverse mapping.
OGRFieldTypes = {
0: OFTInteger,
1: OFTIntegerList,
2: OFTReal,
3: OFTRealList,
4: OFTString,
5: OFTStringList,
6: OFTWideString,
7: OFTWideStringList,
8: OFTBinary,
9: OFTDate,
10: OFTTime,
11: OFTDateTime,
12: OFTInteger64,
13: OFTInteger64List,
}
ROGRFieldTypes = {cls: num for num, cls in OGRFieldTypes.items()}

View File

@@ -0,0 +1,743 @@
"""
The OGRGeometry is a wrapper for using the OGR Geometry class
(see https://gdal.org/api/ogrgeometry_cpp.html#_CPPv411OGRGeometry).
OGRGeometry may be instantiated when reading geometries from OGR Data Sources
(e.g. SHP files), or when given OGC WKT (a string).
While the 'full' API is not present yet, the API is "pythonic" unlike
the traditional and "next-generation" OGR Python bindings. One major
advantage OGR Geometries have over their GEOS counterparts is support
for spatial reference systems and their transformation.
Example:
>>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference
>>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)'
>>> pnt = OGRGeometry(wkt1)
>>> print(pnt)
POINT (-90 30)
>>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84'))
>>> mpnt.add(wkt1)
>>> mpnt.add(wkt1)
>>> print(mpnt)
MULTIPOINT (-90 30,-90 30)
>>> print(mpnt.srs.name)
WGS 84
>>> print(mpnt.srs.proj)
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
>>> mpnt.transform(SpatialReference('NAD27'))
>>> print(mpnt.proj)
+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs
>>> print(mpnt)
MULTIPOINT (-89.99993037860248 29.99979788655764,-89.99993037860248 29.99979788655764)
The OGRGeomType class is to make it easy to specify an OGR geometry type:
>>> from django.contrib.gis.gdal import OGRGeomType
>>> gt1 = OGRGeomType(3) # Using an integer for the type
>>> gt2 = OGRGeomType('Polygon') # Using a string
>>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive
>>> print(gt1 == 3, gt1 == 'Polygon') # Equivalence works w/non-OGRGeomType objects
True True
"""
import sys
from binascii import b2a_hex
from ctypes import byref, c_char_p, c_double, c_ubyte, c_void_p, string_at
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope
from django.contrib.gis.gdal.error import GDALException, SRSException
from django.contrib.gis.gdal.geomtype import OGRGeomType
from django.contrib.gis.gdal.prototypes import geom as capi
from django.contrib.gis.gdal.prototypes import srs as srs_api
from django.contrib.gis.gdal.srs import CoordTransform, SpatialReference
from django.contrib.gis.geometry import hex_regex, json_regex, wkt_regex
from django.utils.encoding import force_bytes
# For more information, see the OGR C API source code:
# https://gdal.org/api/vector_c_api.html
#
# The OGR_G_* routines are relevant here.
class OGRGeometry(GDALBase):
"""Encapsulate an OGR geometry."""
destructor = capi.destroy_geom
def __init__(self, geom_input, srs=None):
"""Initialize Geometry on either WKT or an OGR pointer as input."""
str_instance = isinstance(geom_input, str)
# If HEX, unpack input to a binary buffer.
if str_instance and hex_regex.match(geom_input):
geom_input = memoryview(bytes.fromhex(geom_input))
str_instance = False
# Constructing the geometry,
if str_instance:
wkt_m = wkt_regex.match(geom_input)
json_m = json_regex.match(geom_input)
if wkt_m:
if wkt_m["srid"]:
# If there's EWKT, set the SRS w/value of the SRID.
srs = int(wkt_m["srid"])
if wkt_m["type"].upper() == "LINEARRING":
# OGR_G_CreateFromWkt doesn't work with LINEARRING WKT.
# See https://trac.osgeo.org/gdal/ticket/1992.
g = capi.create_geom(OGRGeomType(wkt_m["type"]).num)
capi.import_wkt(g, byref(c_char_p(wkt_m["wkt"].encode())))
else:
g = capi.from_wkt(
byref(c_char_p(wkt_m["wkt"].encode())), None, byref(c_void_p())
)
elif json_m:
g = self._from_json(geom_input.encode())
else:
# Seeing if the input is a valid short-hand string
# (e.g., 'Point', 'POLYGON').
OGRGeomType(geom_input)
g = capi.create_geom(OGRGeomType(geom_input).num)
elif isinstance(geom_input, memoryview):
# WKB was passed in
g = self._from_wkb(geom_input)
elif isinstance(geom_input, OGRGeomType):
# OGRGeomType was passed in, an empty geometry will be created.
g = capi.create_geom(geom_input.num)
elif isinstance(geom_input, self.ptr_type):
# OGR pointer (c_void_p) was the input.
g = geom_input
else:
raise GDALException(
"Invalid input type for OGR Geometry construction: %s"
% type(geom_input)
)
# Now checking the Geometry pointer before finishing initialization
# by setting the pointer for the object.
if not g:
raise GDALException(
"Cannot create OGR Geometry from input: %s" % geom_input
)
self.ptr = g
# Assigning the SpatialReference object to the geometry, if valid.
if srs:
self.srs = srs
# Setting the class depending upon the OGR Geometry Type
self.__class__ = GEO_CLASSES[self.geom_type.num]
# Pickle routines
def __getstate__(self):
srs = self.srs
if srs:
srs = srs.wkt
else:
srs = None
return bytes(self.wkb), srs
def __setstate__(self, state):
wkb, srs = state
ptr = capi.from_wkb(wkb, None, byref(c_void_p()), len(wkb))
if not ptr:
raise GDALException("Invalid OGRGeometry loaded from pickled state.")
self.ptr = ptr
self.srs = srs
@classmethod
def _from_wkb(cls, geom_input):
return capi.from_wkb(
bytes(geom_input), None, byref(c_void_p()), len(geom_input)
)
@staticmethod
def _from_json(geom_input):
return capi.from_json(geom_input)
@classmethod
def from_bbox(cls, bbox):
"Construct a Polygon from a bounding box (4-tuple)."
x0, y0, x1, y1 = bbox
return OGRGeometry(
"POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))"
% (x0, y0, x0, y1, x1, y1, x1, y0, x0, y0)
)
@staticmethod
def from_json(geom_input):
return OGRGeometry(OGRGeometry._from_json(force_bytes(geom_input)))
@classmethod
def from_gml(cls, gml_string):
return cls(capi.from_gml(force_bytes(gml_string)))
# ### Geometry set-like operations ###
# g = g1 | g2
def __or__(self, other):
"Return the union of the two geometries."
return self.union(other)
# g = g1 & g2
def __and__(self, other):
"Return the intersection of this Geometry and the other."
return self.intersection(other)
# g = g1 - g2
def __sub__(self, other):
"Return the difference this Geometry and the other."
return self.difference(other)
# g = g1 ^ g2
def __xor__(self, other):
"Return the symmetric difference of this Geometry and the other."
return self.sym_difference(other)
def __eq__(self, other):
"Is this Geometry equal to the other?"
return isinstance(other, OGRGeometry) and self.equals(other)
def __str__(self):
"WKT is used for the string representation."
return self.wkt
# #### Geometry Properties ####
@property
def dimension(self):
"Return 0 for points, 1 for lines, and 2 for surfaces."
return capi.get_dims(self.ptr)
def _get_coord_dim(self):
"Return the coordinate dimension of the Geometry."
return capi.get_coord_dim(self.ptr)
def _set_coord_dim(self, dim):
"Set the coordinate dimension of this Geometry."
if dim not in (2, 3):
raise ValueError("Geometry dimension must be either 2 or 3")
capi.set_coord_dim(self.ptr, dim)
coord_dim = property(_get_coord_dim, _set_coord_dim)
@property
def geom_count(self):
"Return the number of elements in this Geometry."
return capi.get_geom_count(self.ptr)
@property
def point_count(self):
"Return the number of Points in this Geometry."
return capi.get_point_count(self.ptr)
@property
def num_points(self):
"Alias for `point_count` (same name method in GEOS API.)"
return self.point_count
@property
def num_coords(self):
"Alias for `point_count`."
return self.point_count
@property
def geom_type(self):
"Return the Type for this Geometry."
return OGRGeomType(capi.get_geom_type(self.ptr))
@property
def geom_name(self):
"Return the Name of this Geometry."
return capi.get_geom_name(self.ptr)
@property
def area(self):
"Return the area for a LinearRing, Polygon, or MultiPolygon; 0 otherwise."
return capi.get_area(self.ptr)
@property
def envelope(self):
"Return the envelope for this Geometry."
# TODO: Fix Envelope() for Point geometries.
return Envelope(capi.get_envelope(self.ptr, byref(OGREnvelope())))
@property
def empty(self):
return capi.is_empty(self.ptr)
@property
def extent(self):
"Return the envelope as a 4-tuple, instead of as an Envelope object."
return self.envelope.tuple
# #### SpatialReference-related Properties ####
# The SRS property
def _get_srs(self):
"Return the Spatial Reference for this Geometry."
try:
srs_ptr = capi.get_geom_srs(self.ptr)
return SpatialReference(srs_api.clone_srs(srs_ptr))
except SRSException:
return None
def _set_srs(self, srs):
"Set the SpatialReference for this geometry."
# Do not have to clone the `SpatialReference` object pointer because
# when it is assigned to this `OGRGeometry` it's internal OGR
# reference count is incremented, and will likewise be released
# (decremented) when this geometry's destructor is called.
if isinstance(srs, SpatialReference):
srs_ptr = srs.ptr
elif isinstance(srs, (int, str)):
sr = SpatialReference(srs)
srs_ptr = sr.ptr
elif srs is None:
srs_ptr = None
else:
raise TypeError(
"Cannot assign spatial reference with object of type: %s" % type(srs)
)
capi.assign_srs(self.ptr, srs_ptr)
srs = property(_get_srs, _set_srs)
# The SRID property
def _get_srid(self):
srs = self.srs
if srs:
return srs.srid
return None
def _set_srid(self, srid):
if isinstance(srid, int) or srid is None:
self.srs = srid
else:
raise TypeError("SRID must be set with an integer.")
srid = property(_get_srid, _set_srid)
# #### Output Methods ####
def _geos_ptr(self):
from django.contrib.gis.geos import GEOSGeometry
return GEOSGeometry._from_wkb(self.wkb)
@property
def geos(self):
"Return a GEOSGeometry object from this OGRGeometry."
from django.contrib.gis.geos import GEOSGeometry
return GEOSGeometry(self._geos_ptr(), self.srid)
@property
def gml(self):
"Return the GML representation of the Geometry."
return capi.to_gml(self.ptr)
@property
def hex(self):
"Return the hexadecimal representation of the WKB (a string)."
return b2a_hex(self.wkb).upper()
@property
def json(self):
"""
Return the GeoJSON representation of this Geometry.
"""
return capi.to_json(self.ptr)
geojson = json
@property
def kml(self):
"Return the KML representation of the Geometry."
return capi.to_kml(self.ptr, None)
@property
def wkb_size(self):
"Return the size of the WKB buffer."
return capi.get_wkbsize(self.ptr)
@property
def wkb(self):
"Return the WKB representation of the Geometry."
if sys.byteorder == "little":
byteorder = 1 # wkbNDR (from ogr_core.h)
else:
byteorder = 0 # wkbXDR
sz = self.wkb_size
# Creating the unsigned character buffer, and passing it in by reference.
buf = (c_ubyte * sz)()
capi.to_wkb(self.ptr, byteorder, byref(buf))
# Returning a buffer of the string at the pointer.
return memoryview(string_at(buf, sz))
@property
def wkt(self):
"Return the WKT representation of the Geometry."
return capi.to_wkt(self.ptr, byref(c_char_p()))
@property
def ewkt(self):
"Return the EWKT representation of the Geometry."
srs = self.srs
if srs and srs.srid:
return "SRID=%s;%s" % (srs.srid, self.wkt)
else:
return self.wkt
# #### Geometry Methods ####
def clone(self):
"Clone this OGR Geometry."
return OGRGeometry(capi.clone_geom(self.ptr), self.srs)
def close_rings(self):
"""
If there are any rings within this geometry that have not been
closed, this routine will do so by adding the starting point at the
end.
"""
# Closing the open rings.
capi.geom_close_rings(self.ptr)
def transform(self, coord_trans, clone=False):
"""
Transform this geometry to a different spatial reference system.
May take a CoordTransform object, a SpatialReference object, string
WKT or PROJ, and/or an integer SRID. By default, return nothing
and transform the geometry in-place. However, if the `clone` keyword is
set, return a transformed clone of this geometry.
"""
if clone:
klone = self.clone()
klone.transform(coord_trans)
return klone
# Depending on the input type, use the appropriate OGR routine
# to perform the transformation.
if isinstance(coord_trans, CoordTransform):
capi.geom_transform(self.ptr, coord_trans.ptr)
elif isinstance(coord_trans, SpatialReference):
capi.geom_transform_to(self.ptr, coord_trans.ptr)
elif isinstance(coord_trans, (int, str)):
sr = SpatialReference(coord_trans)
capi.geom_transform_to(self.ptr, sr.ptr)
else:
raise TypeError(
"Transform only accepts CoordTransform, "
"SpatialReference, string, and integer objects."
)
# #### Topology Methods ####
def _topology(self, func, other):
"""A generalized function for topology operations, takes a GDAL function and
the other geometry to perform the operation on."""
if not isinstance(other, OGRGeometry):
raise TypeError(
"Must use another OGRGeometry object for topology operations!"
)
# Returning the output of the given function with the other geometry's
# pointer.
return func(self.ptr, other.ptr)
def intersects(self, other):
"Return True if this geometry intersects with the other."
return self._topology(capi.ogr_intersects, other)
def equals(self, other):
"Return True if this geometry is equivalent to the other."
return self._topology(capi.ogr_equals, other)
def disjoint(self, other):
"Return True if this geometry and the other are spatially disjoint."
return self._topology(capi.ogr_disjoint, other)
def touches(self, other):
"Return True if this geometry touches the other."
return self._topology(capi.ogr_touches, other)
def crosses(self, other):
"Return True if this geometry crosses the other."
return self._topology(capi.ogr_crosses, other)
def within(self, other):
"Return True if this geometry is within the other."
return self._topology(capi.ogr_within, other)
def contains(self, other):
"Return True if this geometry contains the other."
return self._topology(capi.ogr_contains, other)
def overlaps(self, other):
"Return True if this geometry overlaps the other."
return self._topology(capi.ogr_overlaps, other)
# #### Geometry-generation Methods ####
def _geomgen(self, gen_func, other=None):
"A helper routine for the OGR routines that generate geometries."
if isinstance(other, OGRGeometry):
return OGRGeometry(gen_func(self.ptr, other.ptr), self.srs)
else:
return OGRGeometry(gen_func(self.ptr), self.srs)
@property
def boundary(self):
"Return the boundary of this geometry."
return self._geomgen(capi.get_boundary)
@property
def convex_hull(self):
"""
Return the smallest convex Polygon that contains all the points in
this Geometry.
"""
return self._geomgen(capi.geom_convex_hull)
def difference(self, other):
"""
Return a new geometry consisting of the region which is the difference
of this geometry and the other.
"""
return self._geomgen(capi.geom_diff, other)
def intersection(self, other):
"""
Return a new geometry consisting of the region of intersection of this
geometry and the other.
"""
return self._geomgen(capi.geom_intersection, other)
def sym_difference(self, other):
"""
Return a new geometry which is the symmetric difference of this
geometry and the other.
"""
return self._geomgen(capi.geom_sym_diff, other)
def union(self, other):
"""
Return a new geometry consisting of the region which is the union of
this geometry and the other.
"""
return self._geomgen(capi.geom_union, other)
# The subclasses for OGR Geometry.
class Point(OGRGeometry):
def _geos_ptr(self):
from django.contrib.gis import geos
return geos.Point._create_empty() if self.empty else super()._geos_ptr()
@classmethod
def _create_empty(cls):
return capi.create_geom(OGRGeomType("point").num)
@property
def x(self):
"Return the X coordinate for this Point."
return capi.getx(self.ptr, 0)
@property
def y(self):
"Return the Y coordinate for this Point."
return capi.gety(self.ptr, 0)
@property
def z(self):
"Return the Z coordinate for this Point."
if self.coord_dim == 3:
return capi.getz(self.ptr, 0)
@property
def tuple(self):
"Return the tuple of this point."
if self.coord_dim == 2:
return (self.x, self.y)
elif self.coord_dim == 3:
return (self.x, self.y, self.z)
coords = tuple
class LineString(OGRGeometry):
def __getitem__(self, index):
"Return the Point at the given index."
if 0 <= index < self.point_count:
x, y, z = c_double(), c_double(), c_double()
capi.get_point(self.ptr, index, byref(x), byref(y), byref(z))
dim = self.coord_dim
if dim == 1:
return (x.value,)
elif dim == 2:
return (x.value, y.value)
elif dim == 3:
return (x.value, y.value, z.value)
else:
raise IndexError(
"Index out of range when accessing points of a line string: %s." % index
)
def __len__(self):
"Return the number of points in the LineString."
return self.point_count
@property
def tuple(self):
"Return the tuple representation of this LineString."
return tuple(self[i] for i in range(len(self)))
coords = tuple
def _listarr(self, func):
"""
Internal routine that returns a sequence (list) corresponding with
the given function.
"""
return [func(self.ptr, i) for i in range(len(self))]
@property
def x(self):
"Return the X coordinates in a list."
return self._listarr(capi.getx)
@property
def y(self):
"Return the Y coordinates in a list."
return self._listarr(capi.gety)
@property
def z(self):
"Return the Z coordinates in a list."
if self.coord_dim == 3:
return self._listarr(capi.getz)
# LinearRings are used in Polygons.
class LinearRing(LineString):
pass
class Polygon(OGRGeometry):
def __len__(self):
"Return the number of interior rings in this Polygon."
return self.geom_count
def __getitem__(self, index):
"Get the ring at the specified index."
if 0 <= index < self.geom_count:
return OGRGeometry(
capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs
)
else:
raise IndexError(
"Index out of range when accessing rings of a polygon: %s." % index
)
# Polygon Properties
@property
def shell(self):
"Return the shell of this Polygon."
return self[0] # First ring is the shell
exterior_ring = shell
@property
def tuple(self):
"Return a tuple of LinearRing coordinate tuples."
return tuple(self[i].tuple for i in range(self.geom_count))
coords = tuple
@property
def point_count(self):
"Return the number of Points in this Polygon."
# Summing up the number of points in each ring of the Polygon.
return sum(self[i].point_count for i in range(self.geom_count))
@property
def centroid(self):
"Return the centroid (a Point) of this Polygon."
# The centroid is a Point, create a geometry for this.
p = OGRGeometry(OGRGeomType("Point"))
capi.get_centroid(self.ptr, p.ptr)
return p
# Geometry Collection base class.
class GeometryCollection(OGRGeometry):
"The Geometry Collection class."
def __getitem__(self, index):
"Get the Geometry at the specified index."
if 0 <= index < self.geom_count:
return OGRGeometry(
capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs
)
else:
raise IndexError(
"Index out of range when accessing geometry in a collection: %s."
% index
)
def __len__(self):
"Return the number of geometries in this Geometry Collection."
return self.geom_count
def add(self, geom):
"Add the geometry to this Geometry Collection."
if isinstance(geom, OGRGeometry):
if isinstance(geom, self.__class__):
for g in geom:
capi.add_geom(self.ptr, g.ptr)
else:
capi.add_geom(self.ptr, geom.ptr)
elif isinstance(geom, str):
tmp = OGRGeometry(geom)
capi.add_geom(self.ptr, tmp.ptr)
else:
raise GDALException("Must add an OGRGeometry.")
@property
def point_count(self):
"Return the number of Points in this Geometry Collection."
# Summing up the number of points in each geometry in this collection
return sum(self[i].point_count for i in range(self.geom_count))
@property
def tuple(self):
"Return a tuple representation of this Geometry Collection."
return tuple(self[i].tuple for i in range(self.geom_count))
coords = tuple
# Multiple Geometry types.
class MultiPoint(GeometryCollection):
pass
class MultiLineString(GeometryCollection):
pass
class MultiPolygon(GeometryCollection):
pass
# Class mapping dictionary (using the OGRwkbGeometryType as the key)
GEO_CLASSES = {
1: Point,
2: LineString,
3: Polygon,
4: MultiPoint,
5: MultiLineString,
6: MultiPolygon,
7: GeometryCollection,
101: LinearRing,
1 + OGRGeomType.wkb25bit: Point,
2 + OGRGeomType.wkb25bit: LineString,
3 + OGRGeomType.wkb25bit: Polygon,
4 + OGRGeomType.wkb25bit: MultiPoint,
5 + OGRGeomType.wkb25bit: MultiLineString,
6 + OGRGeomType.wkb25bit: MultiPolygon,
7 + OGRGeomType.wkb25bit: GeometryCollection,
}

View File

@@ -0,0 +1,95 @@
from django.contrib.gis.gdal.error import GDALException
class OGRGeomType:
"Encapsulate OGR Geometry Types."
wkb25bit = -2147483648
# Dictionary of acceptable OGRwkbGeometryType s and their string names.
_types = {
0: "Unknown",
1: "Point",
2: "LineString",
3: "Polygon",
4: "MultiPoint",
5: "MultiLineString",
6: "MultiPolygon",
7: "GeometryCollection",
100: "None",
101: "LinearRing",
102: "PointZ",
1 + wkb25bit: "Point25D",
2 + wkb25bit: "LineString25D",
3 + wkb25bit: "Polygon25D",
4 + wkb25bit: "MultiPoint25D",
5 + wkb25bit: "MultiLineString25D",
6 + wkb25bit: "MultiPolygon25D",
7 + wkb25bit: "GeometryCollection25D",
}
# Reverse type dictionary, keyed by lowercase of the name.
_str_types = {v.lower(): k for k, v in _types.items()}
def __init__(self, type_input):
"Figure out the correct OGR Type based upon the input."
if isinstance(type_input, OGRGeomType):
num = type_input.num
elif isinstance(type_input, str):
type_input = type_input.lower()
if type_input == "geometry":
type_input = "unknown"
num = self._str_types.get(type_input)
if num is None:
raise GDALException('Invalid OGR String Type "%s"' % type_input)
elif isinstance(type_input, int):
if type_input not in self._types:
raise GDALException("Invalid OGR Integer Type: %d" % type_input)
num = type_input
else:
raise TypeError("Invalid OGR input type given.")
# Setting the OGR geometry type number.
self.num = num
def __str__(self):
"Return the value of the name property."
return self.name
def __eq__(self, other):
"""
Do 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.num == other.num
elif isinstance(other, str):
return self.name.lower() == other.lower()
elif isinstance(other, int):
return self.num == other
else:
return False
@property
def name(self):
"Return a short-hand string form of the OGR Geometry type."
return self._types[self.num]
@property
def django(self):
"Return the Django GeometryField for this OGR Type."
s = self.name.replace("25D", "")
if s in ("LinearRing", "None"):
return None
elif s == "Unknown":
s = "Geometry"
elif s == "PointZ":
s = "Point"
return s + "Field"
def to_multi(self):
"""
Transform Point, LineString, Polygon, and their 25D equivalents
to their Multi... counterpart.
"""
if self.name.startswith(("Point", "LineString", "Polygon")):
self.num += 3

View File

@@ -0,0 +1,234 @@
from ctypes import byref, c_double
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.envelope import Envelope, OGREnvelope
from django.contrib.gis.gdal.error import GDALException, SRSException
from django.contrib.gis.gdal.feature import Feature
from django.contrib.gis.gdal.field import OGRFieldTypes
from django.contrib.gis.gdal.geometries import OGRGeometry
from django.contrib.gis.gdal.geomtype import OGRGeomType
from django.contrib.gis.gdal.prototypes import ds as capi
from django.contrib.gis.gdal.prototypes import geom as geom_api
from django.contrib.gis.gdal.prototypes import srs as srs_api
from django.contrib.gis.gdal.srs import SpatialReference
from django.utils.encoding import force_bytes, force_str
# For more information, see the OGR C API source code:
# https://gdal.org/api/vector_c_api.html
#
# The OGR_L_* routines are relevant here.
class Layer(GDALBase):
"""
A class that wraps an OGR Layer, needs to be instantiated from a DataSource
object.
"""
def __init__(self, layer_ptr, ds):
"""
Initialize on an OGR C pointer to the Layer and the `DataSource` object
that owns this layer. The `DataSource` object is required so that a
reference to it is kept with this Layer. This prevents garbage
collection of the `DataSource` while this Layer is still active.
"""
if not layer_ptr:
raise GDALException("Cannot create Layer, invalid pointer given")
self.ptr = layer_ptr
self._ds = ds
self._ldefn = capi.get_layer_defn(self._ptr)
# Does the Layer support random reading?
self._random_read = self.test_capability(b"RandomRead")
def __getitem__(self, index):
"Get the Feature at the specified index."
if isinstance(index, int):
# An integer index was given -- we cannot do a check based on the
# number of features because the beginning and ending feature IDs
# are not guaranteed to be 0 and len(layer)-1, respectively.
if index < 0:
raise IndexError("Negative indices are not allowed on OGR Layers.")
return self._make_feature(index)
elif isinstance(index, slice):
# A slice was given
start, stop, stride = index.indices(self.num_feat)
return [self._make_feature(fid) for fid in range(start, stop, stride)]
else:
raise TypeError(
"Integers and slices may only be used when indexing OGR Layers."
)
def __iter__(self):
"Iterate over each Feature in the Layer."
# ResetReading() must be called before iteration is to begin.
capi.reset_reading(self._ptr)
for i in range(self.num_feat):
yield Feature(capi.get_next_feature(self._ptr), self)
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
def _make_feature(self, feat_id):
"""
Helper routine for __getitem__ that constructs a Feature from the given
Feature ID. If the OGR Layer does not support random-access reading,
then each feature of the layer will be incremented through until the
a Feature is found matching the given feature ID.
"""
if self._random_read:
# If the Layer supports random reading, return.
try:
return Feature(capi.get_feature(self.ptr, feat_id), self)
except GDALException:
pass
else:
# Random access isn't supported, have to increment through
# each feature until the given feature ID is encountered.
for feat in self:
if feat.fid == feat_id:
return feat
# Should have returned a Feature, raise an IndexError.
raise IndexError("Invalid feature id: %s." % feat_id)
# #### Layer properties ####
@property
def extent(self):
"Return the extent (an Envelope) of this layer."
env = OGREnvelope()
capi.get_extent(self.ptr, byref(env), 1)
return Envelope(env)
@property
def name(self):
"Return the name of this layer in the Data Source."
name = capi.get_fd_name(self._ldefn)
return force_str(name, self._ds.encoding, strings_only=True)
@property
def num_feat(self, force=1):
"Return the number of features in the Layer."
return capi.get_feature_count(self.ptr, force)
@property
def num_fields(self):
"Return the number of fields in the Layer."
return capi.get_field_count(self._ldefn)
@property
def geom_type(self):
"Return the geometry type (OGRGeomType) of the Layer."
return OGRGeomType(capi.get_fd_geom_type(self._ldefn))
@property
def srs(self):
"Return the Spatial Reference used in this Layer."
try:
ptr = capi.get_layer_srs(self.ptr)
return SpatialReference(srs_api.clone_srs(ptr))
except SRSException:
return None
@property
def fields(self):
"""
Return a list of string names corresponding to each of the Fields
available in this Layer.
"""
return [
force_str(
capi.get_field_name(capi.get_field_defn(self._ldefn, i)),
self._ds.encoding,
strings_only=True,
)
for i in range(self.num_fields)
]
@property
def field_types(self):
"""
Return a list of the types of fields in this Layer. For example,
return the list [OFTInteger, OFTReal, OFTString] for an OGR layer that
has an integer, a floating-point, and string fields.
"""
return [
OGRFieldTypes[capi.get_field_type(capi.get_field_defn(self._ldefn, i))]
for i in range(self.num_fields)
]
@property
def field_widths(self):
"Return a list of the maximum field widths for the features."
return [
capi.get_field_width(capi.get_field_defn(self._ldefn, i))
for i in range(self.num_fields)
]
@property
def field_precisions(self):
"Return the field precisions for the features."
return [
capi.get_field_precision(capi.get_field_defn(self._ldefn, i))
for i in range(self.num_fields)
]
def _get_spatial_filter(self):
try:
return OGRGeometry(geom_api.clone_geom(capi.get_spatial_filter(self.ptr)))
except GDALException:
return None
def _set_spatial_filter(self, filter):
if isinstance(filter, OGRGeometry):
capi.set_spatial_filter(self.ptr, filter.ptr)
elif isinstance(filter, (tuple, list)):
if not len(filter) == 4:
raise ValueError("Spatial filter list/tuple must have 4 elements.")
# Map c_double onto params -- if a bad type is passed in it
# will be caught here.
xmin, ymin, xmax, ymax = map(c_double, filter)
capi.set_spatial_filter_rect(self.ptr, xmin, ymin, xmax, ymax)
elif filter is None:
capi.set_spatial_filter(self.ptr, None)
else:
raise TypeError(
"Spatial filter must be either an OGRGeometry instance, a 4-tuple, or "
"None."
)
spatial_filter = property(_get_spatial_filter, _set_spatial_filter)
# #### Layer Methods ####
def get_fields(self, field_name):
"""
Return a list containing the given field name for every Feature
in the Layer.
"""
if field_name not in self.fields:
raise GDALException("invalid field name: %s" % field_name)
return [feat.get(field_name) for feat in self]
def get_geoms(self, geos=False):
"""
Return a list containing the OGRGeometry for every Feature in
the Layer.
"""
if geos:
from django.contrib.gis.geos import GEOSGeometry
return [GEOSGeometry(feat.geom.wkb) for feat in self]
else:
return [feat.geom for feat in self]
def test_capability(self, capability):
"""
Return a bool indicating whether the this Layer supports the given
capability (a string). Valid capability strings include:
'RandomRead', 'SequentialWrite', 'RandomWrite', 'FastSpatialFilter',
'FastFeatureCount', 'FastGetExtent', 'CreateField', 'Transactions',
'DeleteFeature', and 'FastSetNextByIndex'.
"""
return bool(capi.test_capability(self.ptr, force_bytes(capability)))

View File

@@ -0,0 +1,138 @@
import logging
import os
import re
from ctypes import CDLL, CFUNCTYPE, c_char_p, c_int
from ctypes.util import find_library
from django.contrib.gis.gdal.error import GDALException
from django.core.exceptions import ImproperlyConfigured
logger = logging.getLogger("django.contrib.gis")
# Custom library path set?
try:
from django.conf import settings
lib_path = settings.GDAL_LIBRARY_PATH
except (AttributeError, ImportError, ImproperlyConfigured, OSError):
lib_path = None
if lib_path:
lib_names = None
elif os.name == "nt":
# Windows NT shared libraries
lib_names = [
"gdal304",
"gdal303",
"gdal302",
"gdal301",
"gdal300",
"gdal204",
"gdal203",
"gdal202",
]
elif os.name == "posix":
# *NIX library names.
lib_names = [
"gdal",
"GDAL",
"gdal3.4.0",
"gdal3.3.0",
"gdal3.2.0",
"gdal3.1.0",
"gdal3.0.0",
"gdal2.4.0",
"gdal2.3.0",
"gdal2.2.0",
]
else:
raise ImproperlyConfigured('GDAL is unsupported on OS "%s".' % os.name)
# Using the ctypes `find_library` utility to find the
# path to the GDAL library from the list of library names.
if lib_names:
for lib_name in lib_names:
lib_path = find_library(lib_name)
if lib_path is not None:
break
if lib_path is None:
raise ImproperlyConfigured(
'Could not find the GDAL library (tried "%s"). Is GDAL installed? '
"If it is, try setting GDAL_LIBRARY_PATH in your settings."
% '", "'.join(lib_names)
)
# This loads the GDAL/OGR C library
lgdal = CDLL(lib_path)
# On Windows, the GDAL binaries have some OSR routines exported with
# STDCALL, while others are not. Thus, the library will also need to
# be loaded up as WinDLL for said OSR functions that require the
# different calling convention.
if os.name == "nt":
from ctypes import WinDLL
lwingdal = WinDLL(lib_path)
def std_call(func):
"""
Return the correct STDCALL function for certain OSR routines on Win32
platforms.
"""
if os.name == "nt":
return lwingdal[func]
else:
return lgdal[func]
# #### Version-information functions. ####
# Return GDAL library version information with the given key.
_version_info = std_call("GDALVersionInfo")
_version_info.argtypes = [c_char_p]
_version_info.restype = c_char_p
def gdal_version():
"Return only the GDAL version number information."
return _version_info(b"RELEASE_NAME")
def gdal_full_version():
"Return the full GDAL version information."
return _version_info(b"")
def gdal_version_info():
ver = gdal_version()
m = re.match(rb"^(?P<major>\d+)\.(?P<minor>\d+)(?:\.(?P<subminor>\d+))?", ver)
if not m:
raise GDALException('Could not parse GDAL version string "%s"' % ver)
major, minor, subminor = m.groups()
return (int(major), int(minor), subminor and int(subminor))
GDAL_VERSION = gdal_version_info()
# Set library error handling so as errors are logged
CPLErrorHandler = CFUNCTYPE(None, c_int, c_int, c_char_p)
def err_handler(error_class, error_number, message):
logger.error("GDAL_ERROR %d: %s", error_number, message)
err_handler = CPLErrorHandler(err_handler)
def function(name, args, restype):
func = std_call(name)
func.argtypes = args
func.restype = restype
return func
set_error_handler = function("CPLSetErrorHandler", [CPLErrorHandler], CPLErrorHandler)
set_error_handler(err_handler)

View File

@@ -0,0 +1,99 @@
"""
This module houses the ctypes function prototypes for OGR DataSource
related data structures. OGR_Dr_*, OGR_DS_*, OGR_L_*, OGR_F_*,
OGR_Fld_* routines are relevant here.
"""
from ctypes import POINTER, c_char_p, c_double, c_int, c_long, c_void_p
from django.contrib.gis.gdal.envelope import OGREnvelope
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.prototypes.generation import (
bool_output,
const_string_output,
double_output,
geom_output,
int64_output,
int_output,
srs_output,
void_output,
voidptr_output,
)
c_int_p = POINTER(c_int) # shortcut type
# Driver Routines
register_all = void_output(lgdal.OGRRegisterAll, [], errcheck=False)
cleanup_all = void_output(lgdal.OGRCleanupAll, [], errcheck=False)
get_driver = voidptr_output(lgdal.OGRGetDriver, [c_int])
get_driver_by_name = voidptr_output(
lgdal.OGRGetDriverByName, [c_char_p], errcheck=False
)
get_driver_count = int_output(lgdal.OGRGetDriverCount, [])
get_driver_name = const_string_output(
lgdal.OGR_Dr_GetName, [c_void_p], decoding="ascii"
)
# DataSource
open_ds = voidptr_output(lgdal.OGROpen, [c_char_p, c_int, POINTER(c_void_p)])
destroy_ds = void_output(lgdal.OGR_DS_Destroy, [c_void_p], errcheck=False)
release_ds = void_output(lgdal.OGRReleaseDataSource, [c_void_p])
get_ds_name = const_string_output(lgdal.OGR_DS_GetName, [c_void_p])
get_layer = voidptr_output(lgdal.OGR_DS_GetLayer, [c_void_p, c_int])
get_layer_by_name = voidptr_output(lgdal.OGR_DS_GetLayerByName, [c_void_p, c_char_p])
get_layer_count = int_output(lgdal.OGR_DS_GetLayerCount, [c_void_p])
# Layer Routines
get_extent = void_output(lgdal.OGR_L_GetExtent, [c_void_p, POINTER(OGREnvelope), c_int])
get_feature = voidptr_output(lgdal.OGR_L_GetFeature, [c_void_p, c_long])
get_feature_count = int_output(lgdal.OGR_L_GetFeatureCount, [c_void_p, c_int])
get_layer_defn = voidptr_output(lgdal.OGR_L_GetLayerDefn, [c_void_p])
get_layer_srs = srs_output(lgdal.OGR_L_GetSpatialRef, [c_void_p])
get_next_feature = voidptr_output(lgdal.OGR_L_GetNextFeature, [c_void_p])
reset_reading = void_output(lgdal.OGR_L_ResetReading, [c_void_p], errcheck=False)
test_capability = int_output(lgdal.OGR_L_TestCapability, [c_void_p, c_char_p])
get_spatial_filter = geom_output(lgdal.OGR_L_GetSpatialFilter, [c_void_p])
set_spatial_filter = void_output(
lgdal.OGR_L_SetSpatialFilter, [c_void_p, c_void_p], errcheck=False
)
set_spatial_filter_rect = void_output(
lgdal.OGR_L_SetSpatialFilterRect,
[c_void_p, c_double, c_double, c_double, c_double],
errcheck=False,
)
# Feature Definition Routines
get_fd_geom_type = int_output(lgdal.OGR_FD_GetGeomType, [c_void_p])
get_fd_name = const_string_output(lgdal.OGR_FD_GetName, [c_void_p])
get_feat_name = const_string_output(lgdal.OGR_FD_GetName, [c_void_p])
get_field_count = int_output(lgdal.OGR_FD_GetFieldCount, [c_void_p])
get_field_defn = voidptr_output(lgdal.OGR_FD_GetFieldDefn, [c_void_p, c_int])
# Feature Routines
clone_feature = voidptr_output(lgdal.OGR_F_Clone, [c_void_p])
destroy_feature = void_output(lgdal.OGR_F_Destroy, [c_void_p], errcheck=False)
feature_equal = int_output(lgdal.OGR_F_Equal, [c_void_p, c_void_p])
get_feat_geom_ref = geom_output(lgdal.OGR_F_GetGeometryRef, [c_void_p])
get_feat_field_count = int_output(lgdal.OGR_F_GetFieldCount, [c_void_p])
get_feat_field_defn = voidptr_output(lgdal.OGR_F_GetFieldDefnRef, [c_void_p, c_int])
get_fid = int_output(lgdal.OGR_F_GetFID, [c_void_p])
get_field_as_datetime = int_output(
lgdal.OGR_F_GetFieldAsDateTime,
[c_void_p, c_int, c_int_p, c_int_p, c_int_p, c_int_p, c_int_p, c_int_p],
)
get_field_as_double = double_output(lgdal.OGR_F_GetFieldAsDouble, [c_void_p, c_int])
get_field_as_integer = int_output(lgdal.OGR_F_GetFieldAsInteger, [c_void_p, c_int])
get_field_as_integer64 = int64_output(
lgdal.OGR_F_GetFieldAsInteger64, [c_void_p, c_int]
)
is_field_set = bool_output(lgdal.OGR_F_IsFieldSetAndNotNull, [c_void_p, c_int])
get_field_as_string = const_string_output(
lgdal.OGR_F_GetFieldAsString, [c_void_p, c_int]
)
get_field_index = int_output(lgdal.OGR_F_GetFieldIndex, [c_void_p, c_char_p])
# Field Routines
get_field_name = const_string_output(lgdal.OGR_Fld_GetNameRef, [c_void_p])
get_field_precision = int_output(lgdal.OGR_Fld_GetPrecision, [c_void_p])
get_field_type = int_output(lgdal.OGR_Fld_GetType, [c_void_p])
get_field_type_name = const_string_output(lgdal.OGR_GetFieldTypeName, [c_int])
get_field_width = int_output(lgdal.OGR_Fld_GetWidth, [c_void_p])

View File

@@ -0,0 +1,141 @@
"""
This module houses the error-checking routines used by the GDAL
ctypes prototypes.
"""
from ctypes import c_void_p, string_at
from django.contrib.gis.gdal.error import GDALException, SRSException, check_err
from django.contrib.gis.gdal.libgdal import lgdal
# Helper routines for retrieving pointers and/or values from
# arguments passed in by reference.
def arg_byref(args, offset=-1):
"Return the pointer argument's by-reference value."
return args[offset]._obj.value
def ptr_byref(args, offset=-1):
"Return the pointer argument passed in by-reference."
return args[offset]._obj
# ### String checking Routines ###
def check_const_string(result, func, cargs, offset=None, cpl=False):
"""
Similar functionality to `check_string`, but does not free the pointer.
"""
if offset:
check_err(result, cpl=cpl)
ptr = ptr_byref(cargs, offset)
return ptr.value
else:
return result
def check_string(result, func, cargs, offset=-1, str_result=False):
"""
Check the string output returned from the given function, and free
the string pointer allocated by OGR. The `str_result` keyword
may be used when the result is the string pointer, otherwise
the OGR error code is assumed. The `offset` keyword may be used
to extract the string pointer passed in by-reference at the given
slice offset in the function arguments.
"""
if str_result:
# For routines that return a string.
ptr = result
if not ptr:
s = None
else:
s = string_at(result)
else:
# Error-code return specified.
check_err(result)
ptr = ptr_byref(cargs, offset)
# Getting the string value
s = ptr.value
# Correctly freeing the allocated memory behind GDAL pointer
# with the VSIFree routine.
if ptr:
lgdal.VSIFree(ptr)
return s
# ### DataSource, Layer error-checking ###
# ### Envelope checking ###
def check_envelope(result, func, cargs, offset=-1):
"Check a function that returns an OGR Envelope by reference."
return ptr_byref(cargs, offset)
# ### Geometry error-checking routines ###
def check_geom(result, func, cargs):
"Check a function that returns a geometry."
# OGR_G_Clone may return an integer, even though the
# restype is set to c_void_p
if isinstance(result, int):
result = c_void_p(result)
if not result:
raise GDALException(
'Invalid geometry pointer returned from "%s".' % func.__name__
)
return result
def check_geom_offset(result, func, cargs, offset=-1):
"Check the geometry at the given offset in the C parameter list."
check_err(result)
geom = ptr_byref(cargs, offset=offset)
return check_geom(geom, func, cargs)
# ### Spatial Reference error-checking routines ###
def check_srs(result, func, cargs):
if isinstance(result, int):
result = c_void_p(result)
if not result:
raise SRSException(
'Invalid spatial reference pointer returned from "%s".' % func.__name__
)
return result
# ### Other error-checking routines ###
def check_arg_errcode(result, func, cargs, cpl=False):
"""
The error code is returned in the last argument, by reference.
Check its value with `check_err` before returning the result.
"""
check_err(arg_byref(cargs), cpl=cpl)
return result
def check_errcode(result, func, cargs, cpl=False):
"""
Check the error code returned (c_int).
"""
check_err(result, cpl=cpl)
def check_pointer(result, func, cargs):
"Make sure the result pointer is valid."
if isinstance(result, int):
result = c_void_p(result)
if result:
return result
else:
raise GDALException('Invalid pointer returned from "%s"' % func.__name__)
def check_str_arg(result, func, cargs):
"""
This is for the OSRGet[Angular|Linear]Units functions, which
require that the returned string pointer not be freed. This
returns both the double and string values.
"""
dbl = result
ptr = cargs[-1]._obj
return dbl, ptr.value.decode()

View File

@@ -0,0 +1,177 @@
"""
This module contains functions that generate ctypes prototypes for the
GDAL routines.
"""
from ctypes import POINTER, c_bool, c_char_p, c_double, c_int, c_int64, c_void_p
from functools import partial
from django.contrib.gis.gdal.prototypes.errcheck import (
check_arg_errcode,
check_const_string,
check_errcode,
check_geom,
check_geom_offset,
check_pointer,
check_srs,
check_str_arg,
check_string,
)
class gdal_char_p(c_char_p):
pass
def bool_output(func, argtypes, errcheck=None):
"""Generate a ctypes function that returns a boolean value."""
func.argtypes = argtypes
func.restype = c_bool
if errcheck:
func.errcheck = errcheck
return func
def double_output(func, argtypes, errcheck=False, strarg=False, cpl=False):
"Generate a ctypes function that returns a double value."
func.argtypes = argtypes
func.restype = c_double
if errcheck:
func.errcheck = partial(check_arg_errcode, cpl=cpl)
if strarg:
func.errcheck = check_str_arg
return func
def geom_output(func, argtypes, offset=None):
"""
Generate a function that returns a Geometry either by reference
or directly (if the return_geom keyword is set to True).
"""
# Setting the argument types
func.argtypes = argtypes
if not offset:
# When a geometry pointer is directly returned.
func.restype = c_void_p
func.errcheck = check_geom
else:
# Error code returned, geometry is returned by-reference.
func.restype = c_int
def geomerrcheck(result, func, cargs):
return check_geom_offset(result, func, cargs, offset)
func.errcheck = geomerrcheck
return func
def int_output(func, argtypes, errcheck=None):
"Generate a ctypes function that returns an integer value."
func.argtypes = argtypes
func.restype = c_int
if errcheck:
func.errcheck = errcheck
return func
def int64_output(func, argtypes):
"Generate a ctypes function that returns a 64-bit integer value."
func.argtypes = argtypes
func.restype = c_int64
return func
def srs_output(func, argtypes):
"""
Generate a ctypes prototype for the given function with
the given C arguments that returns a pointer to an OGR
Spatial Reference System.
"""
func.argtypes = argtypes
func.restype = c_void_p
func.errcheck = check_srs
return func
def const_string_output(func, argtypes, offset=None, decoding=None, cpl=False):
func.argtypes = argtypes
if offset:
func.restype = c_int
else:
func.restype = c_char_p
def _check_const(result, func, cargs):
res = check_const_string(result, func, cargs, offset=offset, cpl=cpl)
if res and decoding:
res = res.decode(decoding)
return res
func.errcheck = _check_const
return func
def string_output(func, argtypes, offset=-1, str_result=False, decoding=None):
"""
Generate a ctypes prototype for the given function with the
given argument types that returns a string from a GDAL pointer.
The `const` flag indicates whether the allocated pointer should
be freed via the GDAL library routine VSIFree -- but only applies
only when `str_result` is True.
"""
func.argtypes = argtypes
if str_result:
# Use subclass of c_char_p so the error checking routine
# can free the memory at the pointer's address.
func.restype = gdal_char_p
else:
# Error code is returned
func.restype = c_int
# Dynamically defining our error-checking function with the
# given offset.
def _check_str(result, func, cargs):
res = check_string(result, func, cargs, offset=offset, str_result=str_result)
if res and decoding:
res = res.decode(decoding)
return res
func.errcheck = _check_str
return func
def void_output(func, argtypes, errcheck=True, cpl=False):
"""
For functions that don't only return an error code that needs to
be examined.
"""
if argtypes:
func.argtypes = argtypes
if errcheck:
# `errcheck` keyword may be set to False for routines that
# return void, rather than a status code.
func.restype = c_int
func.errcheck = partial(check_errcode, cpl=cpl)
else:
func.restype = None
return func
def voidptr_output(func, argtypes, errcheck=True):
"For functions that return c_void_p."
func.argtypes = argtypes
func.restype = c_void_p
if errcheck:
func.errcheck = check_pointer
return func
def chararray_output(func, argtypes, errcheck=True):
"""For functions that return a c_char_p array."""
func.argtypes = argtypes
func.restype = POINTER(c_char_p)
if errcheck:
func.errcheck = check_pointer
return func

View File

@@ -0,0 +1,139 @@
from ctypes import POINTER, c_char_p, c_double, c_int, c_void_p
from django.contrib.gis.gdal.envelope import OGREnvelope
from django.contrib.gis.gdal.libgdal import lgdal
from django.contrib.gis.gdal.prototypes.errcheck import check_envelope
from django.contrib.gis.gdal.prototypes.generation import (
const_string_output,
double_output,
geom_output,
int_output,
srs_output,
string_output,
void_output,
)
# ### Generation routines specific to this module ###
def env_func(f, argtypes):
"For getting OGREnvelopes."
f.argtypes = argtypes
f.restype = None
f.errcheck = check_envelope
return f
def pnt_func(f):
"For accessing point information."
return double_output(f, [c_void_p, c_int])
def topology_func(f):
f.argtypes = [c_void_p, c_void_p]
f.restype = c_int
f.errcheck = lambda result, func, cargs: bool(result)
return f
# ### OGR_G ctypes function prototypes ###
# GeoJSON routines.
from_json = geom_output(lgdal.OGR_G_CreateGeometryFromJson, [c_char_p])
to_json = string_output(
lgdal.OGR_G_ExportToJson, [c_void_p], str_result=True, decoding="ascii"
)
to_kml = string_output(
lgdal.OGR_G_ExportToKML, [c_void_p, c_char_p], str_result=True, decoding="ascii"
)
# GetX, GetY, GetZ all return doubles.
getx = pnt_func(lgdal.OGR_G_GetX)
gety = pnt_func(lgdal.OGR_G_GetY)
getz = pnt_func(lgdal.OGR_G_GetZ)
# Geometry creation routines.
from_wkb = geom_output(
lgdal.OGR_G_CreateFromWkb, [c_char_p, c_void_p, POINTER(c_void_p), c_int], offset=-2
)
from_wkt = geom_output(
lgdal.OGR_G_CreateFromWkt,
[POINTER(c_char_p), c_void_p, POINTER(c_void_p)],
offset=-1,
)
from_gml = geom_output(lgdal.OGR_G_CreateFromGML, [c_char_p])
create_geom = geom_output(lgdal.OGR_G_CreateGeometry, [c_int])
clone_geom = geom_output(lgdal.OGR_G_Clone, [c_void_p])
get_geom_ref = geom_output(lgdal.OGR_G_GetGeometryRef, [c_void_p, c_int])
get_boundary = geom_output(lgdal.OGR_G_GetBoundary, [c_void_p])
geom_convex_hull = geom_output(lgdal.OGR_G_ConvexHull, [c_void_p])
geom_diff = geom_output(lgdal.OGR_G_Difference, [c_void_p, c_void_p])
geom_intersection = geom_output(lgdal.OGR_G_Intersection, [c_void_p, c_void_p])
geom_sym_diff = geom_output(lgdal.OGR_G_SymmetricDifference, [c_void_p, c_void_p])
geom_union = geom_output(lgdal.OGR_G_Union, [c_void_p, c_void_p])
# Geometry modification routines.
add_geom = void_output(lgdal.OGR_G_AddGeometry, [c_void_p, c_void_p])
import_wkt = void_output(lgdal.OGR_G_ImportFromWkt, [c_void_p, POINTER(c_char_p)])
# Destroys a geometry
destroy_geom = void_output(lgdal.OGR_G_DestroyGeometry, [c_void_p], errcheck=False)
# Geometry export routines.
to_wkb = void_output(
lgdal.OGR_G_ExportToWkb, None, errcheck=True
) # special handling for WKB.
to_wkt = string_output(
lgdal.OGR_G_ExportToWkt, [c_void_p, POINTER(c_char_p)], decoding="ascii"
)
to_gml = string_output(
lgdal.OGR_G_ExportToGML, [c_void_p], str_result=True, decoding="ascii"
)
get_wkbsize = int_output(lgdal.OGR_G_WkbSize, [c_void_p])
# Geometry spatial-reference related routines.
assign_srs = void_output(
lgdal.OGR_G_AssignSpatialReference, [c_void_p, c_void_p], errcheck=False
)
get_geom_srs = srs_output(lgdal.OGR_G_GetSpatialReference, [c_void_p])
# Geometry properties
get_area = double_output(lgdal.OGR_G_GetArea, [c_void_p])
get_centroid = void_output(lgdal.OGR_G_Centroid, [c_void_p, c_void_p])
get_dims = int_output(lgdal.OGR_G_GetDimension, [c_void_p])
get_coord_dim = int_output(lgdal.OGR_G_GetCoordinateDimension, [c_void_p])
set_coord_dim = void_output(
lgdal.OGR_G_SetCoordinateDimension, [c_void_p, c_int], errcheck=False
)
is_empty = int_output(
lgdal.OGR_G_IsEmpty, [c_void_p], errcheck=lambda result, func, cargs: bool(result)
)
get_geom_count = int_output(lgdal.OGR_G_GetGeometryCount, [c_void_p])
get_geom_name = const_string_output(
lgdal.OGR_G_GetGeometryName, [c_void_p], decoding="ascii"
)
get_geom_type = int_output(lgdal.OGR_G_GetGeometryType, [c_void_p])
get_point_count = int_output(lgdal.OGR_G_GetPointCount, [c_void_p])
get_point = void_output(
lgdal.OGR_G_GetPoint,
[c_void_p, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double)],
errcheck=False,
)
geom_close_rings = void_output(lgdal.OGR_G_CloseRings, [c_void_p], errcheck=False)
# Topology routines.
ogr_contains = topology_func(lgdal.OGR_G_Contains)
ogr_crosses = topology_func(lgdal.OGR_G_Crosses)
ogr_disjoint = topology_func(lgdal.OGR_G_Disjoint)
ogr_equals = topology_func(lgdal.OGR_G_Equals)
ogr_intersects = topology_func(lgdal.OGR_G_Intersects)
ogr_overlaps = topology_func(lgdal.OGR_G_Overlaps)
ogr_touches = topology_func(lgdal.OGR_G_Touches)
ogr_within = topology_func(lgdal.OGR_G_Within)
# Transformation routines.
geom_transform = void_output(lgdal.OGR_G_Transform, [c_void_p, c_void_p])
geom_transform_to = void_output(lgdal.OGR_G_TransformTo, [c_void_p, c_void_p])
# For retrieving the envelope of the geometry.
get_envelope = env_func(lgdal.OGR_G_GetEnvelope, [c_void_p, POINTER(OGREnvelope)])

View File

@@ -0,0 +1,177 @@
"""
This module houses the ctypes function prototypes for GDAL DataSource (raster)
related data structures.
"""
from ctypes import POINTER, c_bool, c_char_p, c_double, c_int, c_void_p
from functools import partial
from django.contrib.gis.gdal.libgdal import std_call
from django.contrib.gis.gdal.prototypes.generation import (
chararray_output,
const_string_output,
double_output,
int_output,
void_output,
voidptr_output,
)
# For more detail about c function names and definitions see
# https://gdal.org/api/raster_c_api.html
# https://gdal.org/doxygen/gdalwarper_8h.html
# https://gdal.org/api/gdal_utils.html
# Prepare partial functions that use cpl error codes
void_output = partial(void_output, cpl=True)
const_string_output = partial(const_string_output, cpl=True)
double_output = partial(double_output, cpl=True)
# Raster Driver Routines
register_all = void_output(std_call("GDALAllRegister"), [], errcheck=False)
get_driver = voidptr_output(std_call("GDALGetDriver"), [c_int])
get_driver_by_name = voidptr_output(
std_call("GDALGetDriverByName"), [c_char_p], errcheck=False
)
get_driver_count = int_output(std_call("GDALGetDriverCount"), [])
get_driver_description = const_string_output(std_call("GDALGetDescription"), [c_void_p])
# Raster Data Source Routines
create_ds = voidptr_output(
std_call("GDALCreate"), [c_void_p, c_char_p, c_int, c_int, c_int, c_int, c_void_p]
)
open_ds = voidptr_output(std_call("GDALOpen"), [c_char_p, c_int])
close_ds = void_output(std_call("GDALClose"), [c_void_p], errcheck=False)
flush_ds = int_output(std_call("GDALFlushCache"), [c_void_p])
copy_ds = voidptr_output(
std_call("GDALCreateCopy"),
[c_void_p, c_char_p, c_void_p, c_int, POINTER(c_char_p), c_void_p, c_void_p],
)
add_band_ds = void_output(std_call("GDALAddBand"), [c_void_p, c_int])
get_ds_description = const_string_output(std_call("GDALGetDescription"), [c_void_p])
get_ds_driver = voidptr_output(std_call("GDALGetDatasetDriver"), [c_void_p])
get_ds_info = const_string_output(std_call("GDALInfo"), [c_void_p, c_void_p])
get_ds_xsize = int_output(std_call("GDALGetRasterXSize"), [c_void_p])
get_ds_ysize = int_output(std_call("GDALGetRasterYSize"), [c_void_p])
get_ds_raster_count = int_output(std_call("GDALGetRasterCount"), [c_void_p])
get_ds_raster_band = voidptr_output(std_call("GDALGetRasterBand"), [c_void_p, c_int])
get_ds_projection_ref = const_string_output(
std_call("GDALGetProjectionRef"), [c_void_p]
)
set_ds_projection_ref = void_output(std_call("GDALSetProjection"), [c_void_p, c_char_p])
get_ds_geotransform = void_output(
std_call("GDALGetGeoTransform"), [c_void_p, POINTER(c_double * 6)], errcheck=False
)
set_ds_geotransform = void_output(
std_call("GDALSetGeoTransform"), [c_void_p, POINTER(c_double * 6)]
)
get_ds_metadata = chararray_output(
std_call("GDALGetMetadata"), [c_void_p, c_char_p], errcheck=False
)
set_ds_metadata = void_output(
std_call("GDALSetMetadata"), [c_void_p, POINTER(c_char_p), c_char_p]
)
get_ds_metadata_domain_list = chararray_output(
std_call("GDALGetMetadataDomainList"), [c_void_p], errcheck=False
)
get_ds_metadata_item = const_string_output(
std_call("GDALGetMetadataItem"), [c_void_p, c_char_p, c_char_p]
)
set_ds_metadata_item = const_string_output(
std_call("GDALSetMetadataItem"), [c_void_p, c_char_p, c_char_p, c_char_p]
)
free_dsl = void_output(std_call("CSLDestroy"), [POINTER(c_char_p)], errcheck=False)
# Raster Band Routines
band_io = void_output(
std_call("GDALRasterIO"),
[
c_void_p,
c_int,
c_int,
c_int,
c_int,
c_int,
c_void_p,
c_int,
c_int,
c_int,
c_int,
c_int,
],
)
get_band_xsize = int_output(std_call("GDALGetRasterBandXSize"), [c_void_p])
get_band_ysize = int_output(std_call("GDALGetRasterBandYSize"), [c_void_p])
get_band_index = int_output(std_call("GDALGetBandNumber"), [c_void_p])
get_band_description = const_string_output(std_call("GDALGetDescription"), [c_void_p])
get_band_ds = voidptr_output(std_call("GDALGetBandDataset"), [c_void_p])
get_band_datatype = int_output(std_call("GDALGetRasterDataType"), [c_void_p])
get_band_color_interp = int_output(
std_call("GDALGetRasterColorInterpretation"), [c_void_p]
)
get_band_nodata_value = double_output(
std_call("GDALGetRasterNoDataValue"), [c_void_p, POINTER(c_int)]
)
set_band_nodata_value = void_output(
std_call("GDALSetRasterNoDataValue"), [c_void_p, c_double]
)
delete_band_nodata_value = void_output(
std_call("GDALDeleteRasterNoDataValue"), [c_void_p]
)
get_band_statistics = void_output(
std_call("GDALGetRasterStatistics"),
[
c_void_p,
c_int,
c_int,
POINTER(c_double),
POINTER(c_double),
POINTER(c_double),
POINTER(c_double),
c_void_p,
c_void_p,
],
)
compute_band_statistics = void_output(
std_call("GDALComputeRasterStatistics"),
[
c_void_p,
c_int,
POINTER(c_double),
POINTER(c_double),
POINTER(c_double),
POINTER(c_double),
c_void_p,
c_void_p,
],
)
# Reprojection routine
reproject_image = void_output(
std_call("GDALReprojectImage"),
[
c_void_p,
c_char_p,
c_void_p,
c_char_p,
c_int,
c_double,
c_double,
c_void_p,
c_void_p,
c_void_p,
],
)
auto_create_warped_vrt = voidptr_output(
std_call("GDALAutoCreateWarpedVRT"),
[c_void_p, c_char_p, c_char_p, c_int, c_double, c_void_p],
)
# Create VSI gdal raster files from in-memory buffers.
# https://gdal.org/api/cpl.html#cpl-vsi-h
create_vsi_file_from_mem_buffer = voidptr_output(
std_call("VSIFileFromMemBuffer"), [c_char_p, c_void_p, c_int, c_int]
)
get_mem_buffer_from_vsi_file = voidptr_output(
std_call("VSIGetMemFileBuffer"), [c_char_p, POINTER(c_int), c_bool]
)
unlink_vsi_file = int_output(std_call("VSIUnlink"), [c_char_p])

View File

@@ -0,0 +1,109 @@
from ctypes import POINTER, c_char_p, c_int, c_void_p
from django.contrib.gis.gdal.libgdal import GDAL_VERSION, lgdal, std_call
from django.contrib.gis.gdal.prototypes.generation import (
const_string_output,
double_output,
int_output,
srs_output,
string_output,
void_output,
)
# Shortcut generation for routines with known parameters.
def srs_double(f):
"""
Create a function prototype for the OSR routines that take
the OSRSpatialReference object and return a double value.
"""
return double_output(f, [c_void_p, POINTER(c_int)], errcheck=True)
def units_func(f):
"""
Create a ctypes function prototype for OSR units functions, e.g.,
OSRGetAngularUnits, OSRGetLinearUnits.
"""
return double_output(f, [c_void_p, POINTER(c_char_p)], strarg=True)
# Creation & destruction.
clone_srs = srs_output(std_call("OSRClone"), [c_void_p])
new_srs = srs_output(std_call("OSRNewSpatialReference"), [c_char_p])
release_srs = void_output(lgdal.OSRRelease, [c_void_p], errcheck=False)
destroy_srs = void_output(
std_call("OSRDestroySpatialReference"), [c_void_p], errcheck=False
)
srs_validate = void_output(lgdal.OSRValidate, [c_void_p])
if GDAL_VERSION >= (3, 0):
set_axis_strategy = void_output(
lgdal.OSRSetAxisMappingStrategy, [c_void_p, c_int], errcheck=False
)
# Getting the semi_major, semi_minor, and flattening functions.
semi_major = srs_double(lgdal.OSRGetSemiMajor)
semi_minor = srs_double(lgdal.OSRGetSemiMinor)
invflattening = srs_double(lgdal.OSRGetInvFlattening)
# WKT, PROJ, EPSG, XML importation routines.
from_wkt = void_output(lgdal.OSRImportFromWkt, [c_void_p, POINTER(c_char_p)])
from_proj = void_output(lgdal.OSRImportFromProj4, [c_void_p, c_char_p])
from_epsg = void_output(std_call("OSRImportFromEPSG"), [c_void_p, c_int])
from_xml = void_output(lgdal.OSRImportFromXML, [c_void_p, c_char_p])
from_user_input = void_output(std_call("OSRSetFromUserInput"), [c_void_p, c_char_p])
# Morphing to/from ESRI WKT.
morph_to_esri = void_output(lgdal.OSRMorphToESRI, [c_void_p])
morph_from_esri = void_output(lgdal.OSRMorphFromESRI, [c_void_p])
# Identifying the EPSG
identify_epsg = void_output(lgdal.OSRAutoIdentifyEPSG, [c_void_p])
# Getting the angular_units, linear_units functions
linear_units = units_func(lgdal.OSRGetLinearUnits)
angular_units = units_func(lgdal.OSRGetAngularUnits)
# For exporting to WKT, PROJ, "Pretty" WKT, and XML.
to_wkt = string_output(
std_call("OSRExportToWkt"), [c_void_p, POINTER(c_char_p)], decoding="utf-8"
)
to_proj = string_output(
std_call("OSRExportToProj4"), [c_void_p, POINTER(c_char_p)], decoding="ascii"
)
to_pretty_wkt = string_output(
std_call("OSRExportToPrettyWkt"),
[c_void_p, POINTER(c_char_p), c_int],
offset=-2,
decoding="utf-8",
)
to_xml = string_output(
lgdal.OSRExportToXML,
[c_void_p, POINTER(c_char_p), c_char_p],
offset=-2,
decoding="utf-8",
)
# String attribute retrieval routines.
get_attr_value = const_string_output(
std_call("OSRGetAttrValue"), [c_void_p, c_char_p, c_int], decoding="utf-8"
)
get_auth_name = const_string_output(
lgdal.OSRGetAuthorityName, [c_void_p, c_char_p], decoding="ascii"
)
get_auth_code = const_string_output(
lgdal.OSRGetAuthorityCode, [c_void_p, c_char_p], decoding="ascii"
)
# SRS Properties
isgeographic = int_output(lgdal.OSRIsGeographic, [c_void_p])
islocal = int_output(lgdal.OSRIsLocal, [c_void_p])
isprojected = int_output(lgdal.OSRIsProjected, [c_void_p])
# Coordinate transformation
new_ct = srs_output(std_call("OCTNewCoordinateTransformation"), [c_void_p, c_void_p])
destroy_ct = void_output(
std_call("OCTDestroyCoordinateTransformation"), [c_void_p], errcheck=False
)

View File

@@ -0,0 +1,273 @@
from ctypes import byref, c_double, c_int, c_void_p
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.prototypes import raster as capi
from django.contrib.gis.gdal.raster.base import GDALRasterBase
from django.contrib.gis.shortcuts import numpy
from django.utils.encoding import force_str
from .const import (
GDAL_COLOR_TYPES,
GDAL_INTEGER_TYPES,
GDAL_PIXEL_TYPES,
GDAL_TO_CTYPES,
)
class GDALBand(GDALRasterBase):
"""
Wrap a GDAL raster band, needs to be obtained from a GDALRaster object.
"""
def __init__(self, source, index):
self.source = source
self._ptr = capi.get_ds_raster_band(source._ptr, index)
def _flush(self):
"""
Call the flush method on the Band's parent raster and force a refresh
of the statistics attribute when requested the next time.
"""
self.source._flush()
self._stats_refresh = True
@property
def description(self):
"""
Return the description string of the band.
"""
return force_str(capi.get_band_description(self._ptr))
@property
def width(self):
"""
Width (X axis) in pixels of the band.
"""
return capi.get_band_xsize(self._ptr)
@property
def height(self):
"""
Height (Y axis) in pixels of the band.
"""
return capi.get_band_ysize(self._ptr)
@property
def pixel_count(self):
"""
Return the total number of pixels in this band.
"""
return self.width * self.height
_stats_refresh = False
def statistics(self, refresh=False, approximate=False):
"""
Compute statistics on the pixel values of this band.
The return value is a tuple with the following structure:
(minimum, maximum, mean, standard deviation).
If approximate=True, the statistics may be computed based on overviews
or a subset of image tiles.
If refresh=True, the statistics will be computed from the data directly,
and the cache will be updated where applicable.
For empty bands (where all pixel values are nodata), all statistics
values are returned as None.
For raster formats using Persistent Auxiliary Metadata (PAM) services,
the statistics might be cached in an auxiliary file.
"""
# Prepare array with arguments for capi function
smin, smax, smean, sstd = c_double(), c_double(), c_double(), c_double()
stats_args = [
self._ptr,
c_int(approximate),
byref(smin),
byref(smax),
byref(smean),
byref(sstd),
c_void_p(),
c_void_p(),
]
if refresh or self._stats_refresh:
func = capi.compute_band_statistics
else:
# Add additional argument to force computation if there is no
# existing PAM file to take the values from.
force = True
stats_args.insert(2, c_int(force))
func = capi.get_band_statistics
# Computation of statistics fails for empty bands.
try:
func(*stats_args)
result = smin.value, smax.value, smean.value, sstd.value
except GDALException:
result = (None, None, None, None)
self._stats_refresh = False
return result
@property
def min(self):
"""
Return the minimum pixel value for this band.
"""
return self.statistics()[0]
@property
def max(self):
"""
Return the maximum pixel value for this band.
"""
return self.statistics()[1]
@property
def mean(self):
"""
Return the mean of all pixel values of this band.
"""
return self.statistics()[2]
@property
def std(self):
"""
Return the standard deviation of all pixel values of this band.
"""
return self.statistics()[3]
@property
def nodata_value(self):
"""
Return the nodata value for this band, or None if it isn't set.
"""
# Get value and nodata exists flag
nodata_exists = c_int()
value = capi.get_band_nodata_value(self._ptr, nodata_exists)
if not nodata_exists:
value = None
# If the pixeltype is an integer, convert to int
elif self.datatype() in GDAL_INTEGER_TYPES:
value = int(value)
return value
@nodata_value.setter
def nodata_value(self, value):
"""
Set the nodata value for this band.
"""
if value is None:
capi.delete_band_nodata_value(self._ptr)
elif not isinstance(value, (int, float)):
raise ValueError("Nodata value must be numeric or None.")
else:
capi.set_band_nodata_value(self._ptr, value)
self._flush()
def datatype(self, as_string=False):
"""
Return the GDAL Pixel Datatype for this band.
"""
dtype = capi.get_band_datatype(self._ptr)
if as_string:
dtype = GDAL_PIXEL_TYPES[dtype]
return dtype
def color_interp(self, as_string=False):
"""Return the GDAL color interpretation for this band."""
color = capi.get_band_color_interp(self._ptr)
if as_string:
color = GDAL_COLOR_TYPES[color]
return color
def data(self, data=None, offset=None, size=None, shape=None, as_memoryview=False):
"""
Read or writes pixel values for this band. Blocks of data can
be accessed by specifying the width, height and offset of the
desired block. The same specification can be used to update
parts of a raster by providing an array of values.
Allowed input data types are bytes, memoryview, list, tuple, and array.
"""
offset = offset or (0, 0)
size = size or (self.width - offset[0], self.height - offset[1])
shape = shape or size
if any(x <= 0 for x in size):
raise ValueError("Offset too big for this raster.")
if size[0] > self.width or size[1] > self.height:
raise ValueError("Size is larger than raster.")
# Create ctypes type array generator
ctypes_array = GDAL_TO_CTYPES[self.datatype()] * (shape[0] * shape[1])
if data is None:
# Set read mode
access_flag = 0
# Prepare empty ctypes array
data_array = ctypes_array()
else:
# Set write mode
access_flag = 1
# Instantiate ctypes array holding the input data
if isinstance(data, (bytes, memoryview)) or (
numpy and isinstance(data, numpy.ndarray)
):
data_array = ctypes_array.from_buffer_copy(data)
else:
data_array = ctypes_array(*data)
# Access band
capi.band_io(
self._ptr,
access_flag,
offset[0],
offset[1],
size[0],
size[1],
byref(data_array),
shape[0],
shape[1],
self.datatype(),
0,
0,
)
# Return data as numpy array if possible, otherwise as list
if data is None:
if as_memoryview:
return memoryview(data_array)
elif numpy:
# reshape() needs a reshape parameter with the height first.
return numpy.frombuffer(
data_array, dtype=numpy.dtype(data_array)
).reshape(tuple(reversed(size)))
else:
return list(data_array)
else:
self._flush()
class BandList(list):
def __init__(self, source):
self.source = source
super().__init__()
def __iter__(self):
for idx in range(1, len(self) + 1):
yield GDALBand(self.source, idx)
def __len__(self):
return capi.get_ds_raster_count(self.source._ptr)
def __getitem__(self, index):
try:
return GDALBand(self.source, index + 1)
except GDALException:
raise GDALException("Unable to get band index %d" % index)

View File

@@ -0,0 +1,77 @@
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.prototypes import raster as capi
class GDALRasterBase(GDALBase):
"""
Attributes that exist on both GDALRaster and GDALBand.
"""
@property
def metadata(self):
"""
Return the metadata for this raster or band. The return value is a
nested dictionary, where the first-level key is the metadata domain and
the second-level is the metadata item names and values for that domain.
"""
# The initial metadata domain list contains the default domain.
# The default is returned if domain name is None.
domain_list = ["DEFAULT"]
# Get additional metadata domains from the raster.
meta_list = capi.get_ds_metadata_domain_list(self._ptr)
if meta_list:
# The number of domains is unknown, so retrieve data until there
# are no more values in the ctypes array.
counter = 0
domain = meta_list[counter]
while domain:
domain_list.append(domain.decode())
counter += 1
domain = meta_list[counter]
# Free domain list array.
capi.free_dsl(meta_list)
# Retrieve metadata values for each domain.
result = {}
for domain in domain_list:
# Get metadata for this domain.
data = capi.get_ds_metadata(
self._ptr,
(None if domain == "DEFAULT" else domain.encode()),
)
if not data:
continue
# The number of metadata items is unknown, so retrieve data until
# there are no more values in the ctypes array.
domain_meta = {}
counter = 0
item = data[counter]
while item:
key, val = item.decode().split("=")
domain_meta[key] = val
counter += 1
item = data[counter]
# The default domain values are returned if domain is None.
result[domain or "DEFAULT"] = domain_meta
return result
@metadata.setter
def metadata(self, value):
"""
Set the metadata. Update only the domains that are contained in the
value dictionary.
"""
# Loop through domains.
for domain, metadata in value.items():
# Set the domain to None for the default, otherwise encode.
domain = None if domain == "DEFAULT" else domain.encode()
# Set each metadata entry separately.
for meta_name, meta_value in metadata.items():
capi.set_ds_metadata_item(
self._ptr,
meta_name.encode(),
meta_value.encode() if meta_value else None,
domain,
)

View File

@@ -0,0 +1,87 @@
"""
GDAL - Constant definitions
"""
from ctypes import c_double, c_float, c_int16, c_int32, c_ubyte, c_uint16, c_uint32
# See https://gdal.org/api/raster_c_api.html#_CPPv412GDALDataType
GDAL_PIXEL_TYPES = {
0: "GDT_Unknown", # Unknown or unspecified type
1: "GDT_Byte", # Eight bit unsigned integer
2: "GDT_UInt16", # Sixteen bit unsigned integer
3: "GDT_Int16", # Sixteen bit signed integer
4: "GDT_UInt32", # Thirty-two bit unsigned integer
5: "GDT_Int32", # Thirty-two bit signed integer
6: "GDT_Float32", # Thirty-two bit floating point
7: "GDT_Float64", # Sixty-four bit floating point
8: "GDT_CInt16", # Complex Int16
9: "GDT_CInt32", # Complex Int32
10: "GDT_CFloat32", # Complex Float32
11: "GDT_CFloat64", # Complex Float64
}
# A list of gdal datatypes that are integers.
GDAL_INTEGER_TYPES = [1, 2, 3, 4, 5]
# Lookup values to convert GDAL pixel type indices into ctypes objects.
# The GDAL band-io works with ctypes arrays to hold data to be written
# or to hold the space for data to be read into. The lookup below helps
# selecting the right ctypes object for a given gdal pixel type.
GDAL_TO_CTYPES = [
None,
c_ubyte,
c_uint16,
c_int16,
c_uint32,
c_int32,
c_float,
c_double,
None,
None,
None,
None,
]
# List of resampling algorithms that can be used to warp a GDALRaster.
GDAL_RESAMPLE_ALGORITHMS = {
"NearestNeighbour": 0,
"Bilinear": 1,
"Cubic": 2,
"CubicSpline": 3,
"Lanczos": 4,
"Average": 5,
"Mode": 6,
}
# See https://gdal.org/api/raster_c_api.html#_CPPv415GDALColorInterp
GDAL_COLOR_TYPES = {
0: "GCI_Undefined", # Undefined, default value, i.e. not known
1: "GCI_GrayIndex", # Grayscale
2: "GCI_PaletteIndex", # Paletted
3: "GCI_RedBand", # Red band of RGBA image
4: "GCI_GreenBand", # Green band of RGBA image
5: "GCI_BlueBand", # Blue band of RGBA image
6: "GCI_AlphaBand", # Alpha (0=transparent, 255=opaque)
7: "GCI_HueBand", # Hue band of HLS image
8: "GCI_SaturationBand", # Saturation band of HLS image
9: "GCI_LightnessBand", # Lightness band of HLS image
10: "GCI_CyanBand", # Cyan band of CMYK image
11: "GCI_MagentaBand", # Magenta band of CMYK image
12: "GCI_YellowBand", # Yellow band of CMYK image
13: "GCI_BlackBand", # Black band of CMLY image
14: "GCI_YCbCr_YBand", # Y Luminance
15: "GCI_YCbCr_CbBand", # Cb Chroma
16: "GCI_YCbCr_CrBand", # Cr Chroma, also GCI_Max
}
# GDAL virtual filesystems prefix.
VSI_FILESYSTEM_PREFIX = "/vsi"
# Fixed base path for buffer-based GDAL in-memory files.
VSI_MEM_FILESYSTEM_BASE_PATH = "/vsimem/"
# Should the memory file system take ownership of the buffer, freeing it when
# the file is deleted? (No, GDALRaster.__del__() will delete the buffer.)
VSI_TAKE_BUFFER_OWNERSHIP = False
# Should a VSI file be removed when retrieving its buffer?
VSI_DELETE_BUFFER_ON_READ = False

View File

@@ -0,0 +1,539 @@
import json
import os
import sys
import uuid
from ctypes import (
addressof,
byref,
c_buffer,
c_char_p,
c_double,
c_int,
c_void_p,
string_at,
)
from django.contrib.gis.gdal.driver import Driver
from django.contrib.gis.gdal.error import GDALException
from django.contrib.gis.gdal.prototypes import raster as capi
from django.contrib.gis.gdal.raster.band import BandList
from django.contrib.gis.gdal.raster.base import GDALRasterBase
from django.contrib.gis.gdal.raster.const import (
GDAL_RESAMPLE_ALGORITHMS,
VSI_DELETE_BUFFER_ON_READ,
VSI_FILESYSTEM_PREFIX,
VSI_MEM_FILESYSTEM_BASE_PATH,
VSI_TAKE_BUFFER_OWNERSHIP,
)
from django.contrib.gis.gdal.srs import SpatialReference, SRSException
from django.contrib.gis.geometry import json_regex
from django.utils.encoding import force_bytes, force_str
from django.utils.functional import cached_property
class TransformPoint(list):
indices = {
"origin": (0, 3),
"scale": (1, 5),
"skew": (2, 4),
}
def __init__(self, raster, prop):
x = raster.geotransform[self.indices[prop][0]]
y = raster.geotransform[self.indices[prop][1]]
super().__init__([x, y])
self._raster = raster
self._prop = prop
@property
def x(self):
return self[0]
@x.setter
def x(self, value):
gtf = self._raster.geotransform
gtf[self.indices[self._prop][0]] = value
self._raster.geotransform = gtf
@property
def y(self):
return self[1]
@y.setter
def y(self, value):
gtf = self._raster.geotransform
gtf[self.indices[self._prop][1]] = value
self._raster.geotransform = gtf
class GDALRaster(GDALRasterBase):
"""
Wrap a raster GDAL Data Source object.
"""
destructor = capi.close_ds
def __init__(self, ds_input, write=False):
self._write = 1 if write else 0
Driver.ensure_registered()
# Preprocess json inputs. This converts json strings to dictionaries,
# which are parsed below the same way as direct dictionary inputs.
if isinstance(ds_input, str) and json_regex.match(ds_input):
ds_input = json.loads(ds_input)
# If input is a valid file path, try setting file as source.
if isinstance(ds_input, str):
if not ds_input.startswith(VSI_FILESYSTEM_PREFIX) and not os.path.exists(
ds_input
):
raise GDALException(
'Unable to read raster source input "%s".' % ds_input
)
try:
# GDALOpen will auto-detect the data source type.
self._ptr = capi.open_ds(force_bytes(ds_input), self._write)
except GDALException as err:
raise GDALException(
'Could not open the datasource at "{}" ({}).'.format(ds_input, err)
)
elif isinstance(ds_input, bytes):
# Create a new raster in write mode.
self._write = 1
# Get size of buffer.
size = sys.getsizeof(ds_input)
# Pass data to ctypes, keeping a reference to the ctypes object so
# that the vsimem file remains available until the GDALRaster is
# deleted.
self._ds_input = c_buffer(ds_input)
# Create random name to reference in vsimem filesystem.
vsi_path = os.path.join(VSI_MEM_FILESYSTEM_BASE_PATH, str(uuid.uuid4()))
# Create vsimem file from buffer.
capi.create_vsi_file_from_mem_buffer(
force_bytes(vsi_path),
byref(self._ds_input),
size,
VSI_TAKE_BUFFER_OWNERSHIP,
)
# Open the new vsimem file as a GDALRaster.
try:
self._ptr = capi.open_ds(force_bytes(vsi_path), self._write)
except GDALException:
# Remove the broken file from the VSI filesystem.
capi.unlink_vsi_file(force_bytes(vsi_path))
raise GDALException("Failed creating VSI raster from the input buffer.")
elif isinstance(ds_input, dict):
# A new raster needs to be created in write mode
self._write = 1
# Create driver (in memory by default)
driver = Driver(ds_input.get("driver", "MEM"))
# For out of memory drivers, check filename argument
if driver.name != "MEM" and "name" not in ds_input:
raise GDALException(
'Specify name for creation of raster with driver "{}".'.format(
driver.name
)
)
# Check if width and height where specified
if "width" not in ds_input or "height" not in ds_input:
raise GDALException(
"Specify width and height attributes for JSON or dict input."
)
# Check if srid was specified
if "srid" not in ds_input:
raise GDALException("Specify srid for JSON or dict input.")
# Create null terminated gdal options array.
papsz_options = []
for key, val in ds_input.get("papsz_options", {}).items():
option = "{}={}".format(key, val)
papsz_options.append(option.upper().encode())
papsz_options.append(None)
# Convert papszlist to ctypes array.
papsz_options = (c_char_p * len(papsz_options))(*papsz_options)
# Create GDAL Raster
self._ptr = capi.create_ds(
driver._ptr,
force_bytes(ds_input.get("name", "")),
ds_input["width"],
ds_input["height"],
ds_input.get("nr_of_bands", len(ds_input.get("bands", []))),
ds_input.get("datatype", 6),
byref(papsz_options),
)
# Set band data if provided
for i, band_input in enumerate(ds_input.get("bands", [])):
band = self.bands[i]
if "nodata_value" in band_input:
band.nodata_value = band_input["nodata_value"]
# Instantiate band filled with nodata values if only
# partial input data has been provided.
if band.nodata_value is not None and (
"data" not in band_input
or "size" in band_input
or "shape" in band_input
):
band.data(data=(band.nodata_value,), shape=(1, 1))
# Set band data values from input.
band.data(
data=band_input.get("data"),
size=band_input.get("size"),
shape=band_input.get("shape"),
offset=band_input.get("offset"),
)
# Set SRID
self.srs = ds_input.get("srid")
# Set additional properties if provided
if "origin" in ds_input:
self.origin.x, self.origin.y = ds_input["origin"]
if "scale" in ds_input:
self.scale.x, self.scale.y = ds_input["scale"]
if "skew" in ds_input:
self.skew.x, self.skew.y = ds_input["skew"]
elif isinstance(ds_input, c_void_p):
# Instantiate the object using an existing pointer to a gdal raster.
self._ptr = ds_input
else:
raise GDALException(
'Invalid data source input type: "{}".'.format(type(ds_input))
)
def __del__(self):
if self.is_vsi_based:
# Remove the temporary file from the VSI in-memory filesystem.
capi.unlink_vsi_file(force_bytes(self.name))
super().__del__()
def __str__(self):
return self.name
def __repr__(self):
"""
Short-hand representation because WKB may be very large.
"""
return "<Raster object at %s>" % hex(addressof(self._ptr))
def _flush(self):
"""
Flush all data from memory into the source file if it exists.
The data that needs flushing are geotransforms, coordinate systems,
nodata_values and pixel values. This function will be called
automatically wherever it is needed.
"""
# Raise an Exception if the value is being changed in read mode.
if not self._write:
raise GDALException(
"Raster needs to be opened in write mode to change values."
)
capi.flush_ds(self._ptr)
@property
def vsi_buffer(self):
if not (
self.is_vsi_based and self.name.startswith(VSI_MEM_FILESYSTEM_BASE_PATH)
):
return None
# Prepare an integer that will contain the buffer length.
out_length = c_int()
# Get the data using the vsi file name.
dat = capi.get_mem_buffer_from_vsi_file(
force_bytes(self.name),
byref(out_length),
VSI_DELETE_BUFFER_ON_READ,
)
# Read the full buffer pointer.
return string_at(dat, out_length.value)
@cached_property
def is_vsi_based(self):
return self._ptr and self.name.startswith(VSI_FILESYSTEM_PREFIX)
@property
def name(self):
"""
Return the name of this raster. Corresponds to filename
for file-based rasters.
"""
return force_str(capi.get_ds_description(self._ptr))
@cached_property
def driver(self):
"""
Return the GDAL Driver used for this raster.
"""
ds_driver = capi.get_ds_driver(self._ptr)
return Driver(ds_driver)
@property
def width(self):
"""
Width (X axis) in pixels.
"""
return capi.get_ds_xsize(self._ptr)
@property
def height(self):
"""
Height (Y axis) in pixels.
"""
return capi.get_ds_ysize(self._ptr)
@property
def srs(self):
"""
Return the SpatialReference used in this GDALRaster.
"""
try:
wkt = capi.get_ds_projection_ref(self._ptr)
if not wkt:
return None
return SpatialReference(wkt, srs_type="wkt")
except SRSException:
return None
@srs.setter
def srs(self, value):
"""
Set the spatial reference used in this GDALRaster. The input can be
a SpatialReference or any parameter accepted by the SpatialReference
constructor.
"""
if isinstance(value, SpatialReference):
srs = value
elif isinstance(value, (int, str)):
srs = SpatialReference(value)
else:
raise ValueError("Could not create a SpatialReference from input.")
capi.set_ds_projection_ref(self._ptr, srs.wkt.encode())
self._flush()
@property
def srid(self):
"""
Shortcut to access the srid of this GDALRaster.
"""
return self.srs.srid
@srid.setter
def srid(self, value):
"""
Shortcut to set this GDALRaster's srs from an srid.
"""
self.srs = value
@property
def geotransform(self):
"""
Return the geotransform of the data source.
Return the default geotransform if it does not exist or has not been
set previously. The default is [0.0, 1.0, 0.0, 0.0, 0.0, -1.0].
"""
# Create empty ctypes double array for data
gtf = (c_double * 6)()
capi.get_ds_geotransform(self._ptr, byref(gtf))
return list(gtf)
@geotransform.setter
def geotransform(self, values):
"Set the geotransform for the data source."
if len(values) != 6 or not all(isinstance(x, (int, float)) for x in values):
raise ValueError("Geotransform must consist of 6 numeric values.")
# Create ctypes double array with input and write data
values = (c_double * 6)(*values)
capi.set_ds_geotransform(self._ptr, byref(values))
self._flush()
@property
def origin(self):
"""
Coordinates of the raster origin.
"""
return TransformPoint(self, "origin")
@property
def scale(self):
"""
Pixel scale in units of the raster projection.
"""
return TransformPoint(self, "scale")
@property
def skew(self):
"""
Skew of pixels (rotation parameters).
"""
return TransformPoint(self, "skew")
@property
def extent(self):
"""
Return the extent as a 4-tuple (xmin, ymin, xmax, ymax).
"""
# Calculate boundary values based on scale and size
xval = self.origin.x + self.scale.x * self.width
yval = self.origin.y + self.scale.y * self.height
# Calculate min and max values
xmin = min(xval, self.origin.x)
xmax = max(xval, self.origin.x)
ymin = min(yval, self.origin.y)
ymax = max(yval, self.origin.y)
return xmin, ymin, xmax, ymax
@property
def bands(self):
return BandList(self)
def warp(self, ds_input, resampling="NearestNeighbour", max_error=0.0):
"""
Return a warped GDALRaster with the given input characteristics.
The input is expected to be a dictionary containing the parameters
of the target raster. Allowed values are width, height, SRID, origin,
scale, skew, datatype, driver, and name (filename).
By default, the warp functions keeps all parameters equal to the values
of the original source raster. For the name of the target raster, the
name of the source raster will be used and appended with
_copy. + source_driver_name.
In addition, the resampling algorithm can be specified with the "resampling"
input parameter. The default is NearestNeighbor. For a list of all options
consult the GDAL_RESAMPLE_ALGORITHMS constant.
"""
# Get the parameters defining the geotransform, srid, and size of the raster
ds_input.setdefault("width", self.width)
ds_input.setdefault("height", self.height)
ds_input.setdefault("srid", self.srs.srid)
ds_input.setdefault("origin", self.origin)
ds_input.setdefault("scale", self.scale)
ds_input.setdefault("skew", self.skew)
# Get the driver, name, and datatype of the target raster
ds_input.setdefault("driver", self.driver.name)
if "name" not in ds_input:
ds_input["name"] = self.name + "_copy." + self.driver.name
if "datatype" not in ds_input:
ds_input["datatype"] = self.bands[0].datatype()
# Instantiate raster bands filled with nodata values.
ds_input["bands"] = [{"nodata_value": bnd.nodata_value} for bnd in self.bands]
# Create target raster
target = GDALRaster(ds_input, write=True)
# Select resampling algorithm
algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling]
# Reproject image
capi.reproject_image(
self._ptr,
self.srs.wkt.encode(),
target._ptr,
target.srs.wkt.encode(),
algorithm,
0.0,
max_error,
c_void_p(),
c_void_p(),
c_void_p(),
)
# Make sure all data is written to file
target._flush()
return target
def clone(self, name=None):
"""Return a clone of this GDALRaster."""
if name:
clone_name = name
elif self.driver.name != "MEM":
clone_name = self.name + "_copy." + self.driver.name
else:
clone_name = os.path.join(VSI_MEM_FILESYSTEM_BASE_PATH, str(uuid.uuid4()))
return GDALRaster(
capi.copy_ds(
self.driver._ptr,
force_bytes(clone_name),
self._ptr,
c_int(),
c_char_p(),
c_void_p(),
c_void_p(),
),
write=self._write,
)
def transform(
self, srs, driver=None, name=None, resampling="NearestNeighbour", max_error=0.0
):
"""
Return a copy of this raster reprojected into the given spatial
reference system.
"""
# Convert the resampling algorithm name into an algorithm id
algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling]
if isinstance(srs, SpatialReference):
target_srs = srs
elif isinstance(srs, (int, str)):
target_srs = SpatialReference(srs)
else:
raise TypeError(
"Transform only accepts SpatialReference, string, and integer "
"objects."
)
if target_srs.srid == self.srid and (not driver or driver == self.driver.name):
return self.clone(name)
# Create warped virtual dataset in the target reference system
target = capi.auto_create_warped_vrt(
self._ptr,
self.srs.wkt.encode(),
target_srs.wkt.encode(),
algorithm,
max_error,
c_void_p(),
)
target = GDALRaster(target)
# Construct the target warp dictionary from the virtual raster
data = {
"srid": target_srs.srid,
"width": target.width,
"height": target.height,
"origin": [target.origin.x, target.origin.y],
"scale": [target.scale.x, target.scale.y],
"skew": [target.skew.x, target.skew.y],
}
# Set the driver and filepath if provided
if driver:
data["driver"] = driver
if name:
data["name"] = name
# Warp the raster into new srid
return self.warp(data, resampling=resampling, max_error=max_error)
@property
def info(self):
"""
Return information about this raster in a string format equivalent
to the output of the gdalinfo command line utility.
"""
return capi.get_ds_info(self.ptr, None).decode()

View File

@@ -0,0 +1,360 @@
"""
The Spatial Reference class, represents OGR Spatial Reference objects.
Example:
>>> from django.contrib.gis.gdal import SpatialReference
>>> srs = SpatialReference('WGS84')
>>> print(srs)
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"]]
>>> print(srs.proj)
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
>>> print(srs.ellipsoid)
(6378137.0, 6356752.3142451793, 298.25722356300003)
>>> print(srs.projected, srs.geographic)
False True
>>> srs.import_epsg(32140)
>>> print(srs.name)
NAD83 / Texas South Central
"""
from ctypes import byref, c_char_p, c_int
from enum import IntEnum
from django.contrib.gis.gdal.base import GDALBase
from django.contrib.gis.gdal.error import SRSException
from django.contrib.gis.gdal.libgdal import GDAL_VERSION
from django.contrib.gis.gdal.prototypes import srs as capi
from django.utils.encoding import force_bytes, force_str
class AxisOrder(IntEnum):
TRADITIONAL = 0
AUTHORITY = 1
class SpatialReference(GDALBase):
"""
A wrapper for the OGRSpatialReference object. According to the GDAL web site,
the SpatialReference object "provide[s] services to represent coordinate
systems (projections and datums) and to transform between them."
"""
destructor = capi.release_srs
def __init__(self, srs_input="", srs_type="user", axis_order=None):
"""
Create a GDAL OSR Spatial Reference object from the given input.
The input may be string of OGC Well Known Text (WKT), an integer
EPSG code, a PROJ string, and/or a projection "well known" shorthand
string (one of 'WGS84', 'WGS72', 'NAD27', 'NAD83').
"""
if not isinstance(axis_order, (type(None), AxisOrder)):
raise ValueError(
"SpatialReference.axis_order must be an AxisOrder instance."
)
self.axis_order = axis_order or AxisOrder.TRADITIONAL
if srs_type == "wkt":
self.ptr = capi.new_srs(c_char_p(b""))
self.import_wkt(srs_input)
if self.axis_order == AxisOrder.TRADITIONAL and GDAL_VERSION >= (3, 0):
capi.set_axis_strategy(self.ptr, self.axis_order)
elif self.axis_order != AxisOrder.TRADITIONAL and GDAL_VERSION < (3, 0):
raise ValueError("%s is not supported in GDAL < 3.0." % self.axis_order)
return
elif isinstance(srs_input, str):
try:
# If SRID is a string, e.g., '4326', then make acceptable
# as user input.
srid = int(srs_input)
srs_input = "EPSG:%d" % srid
except ValueError:
pass
elif isinstance(srs_input, int):
# EPSG integer code was input.
srs_type = "epsg"
elif isinstance(srs_input, self.ptr_type):
srs = srs_input
srs_type = "ogr"
else:
raise TypeError('Invalid SRS type "%s"' % srs_type)
if srs_type == "ogr":
# Input is already an SRS pointer.
srs = srs_input
else:
# Creating a new SRS pointer, using the string buffer.
buf = c_char_p(b"")
srs = capi.new_srs(buf)
# If the pointer is NULL, throw an exception.
if not srs:
raise SRSException(
"Could not create spatial reference from: %s" % srs_input
)
else:
self.ptr = srs
if self.axis_order == AxisOrder.TRADITIONAL and GDAL_VERSION >= (3, 0):
capi.set_axis_strategy(self.ptr, self.axis_order)
elif self.axis_order != AxisOrder.TRADITIONAL and GDAL_VERSION < (3, 0):
raise ValueError("%s is not supported in GDAL < 3.0." % self.axis_order)
# Importing from either the user input string or an integer SRID.
if srs_type == "user":
self.import_user_input(srs_input)
elif srs_type == "epsg":
self.import_epsg(srs_input)
def __getitem__(self, target):
"""
Return 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 index of the attribute in the WKT. For example:
>>> wkt = 'GEOGCS["WGS 84", DATUM["WGS_1984, ... AUTHORITY["EPSG","4326"]]'
>>> srs = SpatialReference(wkt) # could also use 'WGS84', or 4326
>>> print(srs['GEOGCS'])
WGS 84
>>> print(srs['DATUM'])
WGS_1984
>>> print(srs['AUTHORITY'])
EPSG
>>> print(srs['AUTHORITY', 1]) # The authority value
4326
>>> print(srs['TOWGS84', 4]) # the fourth value in this wkt
0
>>> # For the units authority, have to use the pipe symbole.
>>> print(srs['UNIT|AUTHORITY'])
EPSG
>>> print(srs['UNIT|AUTHORITY', 1]) # The authority value for the units
9122
"""
if isinstance(target, tuple):
return self.attr_value(*target)
else:
return self.attr_value(target)
def __str__(self):
"Use 'pretty' WKT."
return self.pretty_wkt
# #### SpatialReference Methods ####
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.
"""
if not isinstance(target, str) or not isinstance(index, int):
raise TypeError
return capi.get_attr_value(self.ptr, force_bytes(target), index)
def auth_name(self, target):
"Return the authority name for the given string target node."
return capi.get_auth_name(self.ptr, force_bytes(target))
def auth_code(self, target):
"Return the authority code for the given string target node."
return capi.get_auth_code(self.ptr, force_bytes(target))
def clone(self):
"Return a clone of this SpatialReference object."
return SpatialReference(capi.clone_srs(self.ptr), axis_order=self.axis_order)
def from_esri(self):
"Morph this SpatialReference from ESRI's format to EPSG."
capi.morph_from_esri(self.ptr)
def identify_epsg(self):
"""
This method inspects the WKT of this SpatialReference, and will
add EPSG authority nodes where an EPSG identifier is applicable.
"""
capi.identify_epsg(self.ptr)
def to_esri(self):
"Morph this SpatialReference to ESRI's format."
capi.morph_to_esri(self.ptr)
def validate(self):
"Check to see if the given spatial reference is valid."
capi.srs_validate(self.ptr)
# #### Name & SRID properties ####
@property
def name(self):
"Return 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
@property
def srid(self):
"Return the SRID of top-level authority, or None if undefined."
try:
return int(self.attr_value("AUTHORITY", 1))
except (TypeError, ValueError):
return None
# #### Unit Properties ####
@property
def linear_name(self):
"Return the name of the linear units."
units, name = capi.linear_units(self.ptr, byref(c_char_p()))
return name
@property
def linear_units(self):
"Return the value of the linear units."
units, name = capi.linear_units(self.ptr, byref(c_char_p()))
return units
@property
def angular_name(self):
"Return the name of the angular units."
units, name = capi.angular_units(self.ptr, byref(c_char_p()))
return name
@property
def angular_units(self):
"Return the value of the angular units."
units, name = capi.angular_units(self.ptr, byref(c_char_p()))
return units
@property
def units(self):
"""
Return a 2-tuple of the units value and the units name. Automatically
determine whether to return the linear or angular units.
"""
units, name = None, None
if self.projected or self.local:
units, name = capi.linear_units(self.ptr, byref(c_char_p()))
elif self.geographic:
units, name = capi.angular_units(self.ptr, byref(c_char_p()))
if name is not None:
name = force_str(name)
return (units, name)
# #### Spheroid/Ellipsoid Properties ####
@property
def ellipsoid(self):
"""
Return a tuple of the ellipsoid parameters:
(semimajor axis, semiminor axis, and inverse flattening)
"""
return (self.semi_major, self.semi_minor, self.inverse_flattening)
@property
def semi_major(self):
"Return the Semi Major Axis for this Spatial Reference."
return capi.semi_major(self.ptr, byref(c_int()))
@property
def semi_minor(self):
"Return the Semi Minor Axis for this Spatial Reference."
return capi.semi_minor(self.ptr, byref(c_int()))
@property
def inverse_flattening(self):
"Return the Inverse Flattening for this Spatial Reference."
return capi.invflattening(self.ptr, byref(c_int()))
# #### Boolean Properties ####
@property
def geographic(self):
"""
Return True if this SpatialReference is geographic
(root node is GEOGCS).
"""
return bool(capi.isgeographic(self.ptr))
@property
def local(self):
"Return True if this SpatialReference is local (root node is LOCAL_CS)."
return bool(capi.islocal(self.ptr))
@property
def projected(self):
"""
Return True if this SpatialReference is a projected coordinate system
(root node is PROJCS).
"""
return bool(capi.isprojected(self.ptr))
# #### Import Routines #####
def import_epsg(self, epsg):
"Import the Spatial Reference from the EPSG code (an integer)."
capi.from_epsg(self.ptr, epsg)
def import_proj(self, proj):
"""Import the Spatial Reference from a PROJ string."""
capi.from_proj(self.ptr, proj)
def import_user_input(self, user_input):
"Import the Spatial Reference from the given user input string."
capi.from_user_input(self.ptr, force_bytes(user_input))
def import_wkt(self, wkt):
"Import the Spatial Reference from OGC WKT (string)"
capi.from_wkt(self.ptr, byref(c_char_p(force_bytes(wkt))))
def import_xml(self, xml):
"Import the Spatial Reference from an XML string."
capi.from_xml(self.ptr, xml)
# #### Export Properties ####
@property
def wkt(self):
"Return the WKT representation of this Spatial Reference."
return capi.to_wkt(self.ptr, byref(c_char_p()))
@property
def pretty_wkt(self, simplify=0):
"Return the 'pretty' representation of the WKT."
return capi.to_pretty_wkt(self.ptr, byref(c_char_p()), simplify)
@property
def proj(self):
"""Return the PROJ representation for this Spatial Reference."""
return capi.to_proj(self.ptr, byref(c_char_p()))
@property
def proj4(self):
"Alias for proj()."
return self.proj
@property
def xml(self, dialect=""):
"Return the XML representation of this Spatial Reference."
return capi.to_xml(self.ptr, byref(c_char_p()), force_bytes(dialect))
class CoordTransform(GDALBase):
"The coordinate system transformation object."
destructor = capi.destroy_ct
def __init__(self, source, target):
"Initialize on a source and target SpatialReference objects."
if not isinstance(source, SpatialReference) or not isinstance(
target, SpatialReference
):
raise TypeError("source and target must be of type SpatialReference")
self.ptr = capi.new_ct(source._ptr, target._ptr)
self._srs1_name = source.name
self._srs2_name = target.name
def __str__(self):
return 'Transform from "%s" to "%s"' % (self._srs1_name, self._srs2_name)