17.12

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/envelope.py

@@ -0,0 +1,178 @@
+"""
+ 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://www.gdal.org/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)

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/geometries.py

@@ -0,0 +1,715 @@
+"""
+ The OGRGeometry is a wrapper for using the OGR Geometry class
+ (see https://www.gdal.org/classOGRGeometry.html).  OGRGeometry
+ may be instantiated when reading geometries from OGR Data Sources
+ (e.g. SHP files), or when given OGC WKT (a string).
+
+ While the 'full' API is not present yet, the API is "pythonic" unlike
+ the traditional and "next-generation" OGR Python bindings.  One major
+ advantage OGR Geometries have over their GEOS counterparts is support
+ for spatial reference systems and their transformation.
+
+ Example:
+  >>> from django.contrib.gis.gdal import OGRGeometry, OGRGeomType, SpatialReference
+  >>> wkt1, wkt2 = 'POINT(-90 30)', 'POLYGON((0 0, 5 0, 5 5, 0 5)'
+  >>> pnt = OGRGeometry(wkt1)
+  >>> print(pnt)
+  POINT (-90 30)
+  >>> mpnt = OGRGeometry(OGRGeomType('MultiPoint'), SpatialReference('WGS84'))
+  >>> mpnt.add(wkt1)
+  >>> mpnt.add(wkt1)
+  >>> print(mpnt)
+  MULTIPOINT (-90 30,-90 30)
+  >>> print(mpnt.srs.name)
+  WGS 84
+  >>> print(mpnt.srs.proj)
+  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+  >>> mpnt.transform(SpatialReference('NAD27'))
+  >>> print(mpnt.proj)
+  +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs
+  >>> print(mpnt)
+  MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641)
+
+  The OGRGeomType class is to make it easy to specify an OGR geometry type:
+  >>> from django.contrib.gis.gdal import OGRGeomType
+  >>> gt1 = OGRGeomType(3) # Using an integer for the type
+  >>> gt2 = OGRGeomType('Polygon') # Using a string
+  >>> gt3 = OGRGeomType('POLYGON') # It's case-insensitive
+  >>> print(gt1 == 3, gt1 == 'Polygon') # Equivalence works w/non-OGRGeomType objects
+  True 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.libgdal import GDAL_VERSION
+from django.contrib.gis.gdal.prototypes import geom as capi, 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://www.gdal.org/ogr__api_8h.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.group('srid'):
+                    # If there's EWKT, set the SRS w/value of the SRID.
+                    srs = int(wkt_m.group('srid'))
+                if wkt_m.group('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.group('type')).num)
+                    capi.import_wkt(g, byref(c_char_p(wkt_m.group('wkt').encode())))
+                else:
+                    g = capi.from_wkt(byref(c_char_p(wkt_m.group('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):
+        ptr = capi.from_json(geom_input)
+        if GDAL_VERSION < (2, 0):
+            try:
+                capi.get_geom_srs(ptr)
+            except SRSException:
+                srs = SpatialReference(4326)
+                capi.assign_srs(ptr, srs.ptr)
+        return ptr
+
+    @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.4, 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,
+}

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/prototypes/generation.py

@@ -0,0 +1,158 @@
+"""
+ This module contains functions that generate ctypes prototypes for the
+ GDAL routines.
+"""
+from ctypes import POINTER, 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 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

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/prototypes/raster.py

@@ -0,0 +1,108 @@
+"""
+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 GDAL_VERSION, 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/gdal_8h.html
+# https://gdal.org/gdalwarper_8h.html
+# https://www.gdal.org/gdal__utils_8h.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_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)
+
+if GDAL_VERSION >= (2, 1):
+    get_ds_info = const_string_output(std_call('GDALInfo'), [c_void_p, c_void_p])
+else:
+    get_ds_info = None
+
+# 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])
+if GDAL_VERSION >= (2, 1):
+    delete_band_nodata_value = void_output(std_call('GDALDeleteRasterNoDataValue'), [c_void_p])
+else:
+    delete_band_nodata_value = None
+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/cpl__vsi_8h.html
+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])

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/prototypes/srs.py

@@ -0,0 +1,80 @@
+from ctypes import POINTER, c_char_p, c_int, c_void_p
+
+from django.contrib.gis.gdal.libgdal import 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])
+
+# 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.4, "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'
+)
+
+# Memory leak fixed in GDAL 1.5; still exists in 1.4.
+to_xml = string_output(lgdal.OSRExportToXML, [c_void_p, POINTER(c_char_p), c_char_p], offset=-2, decoding='utf-8')
+
+# String attribute retrival 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)

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/raster/band.py

@@ -0,0 +1,252 @@
+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_text
+
+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_text(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:
+            if not capi.delete_band_nodata_value:
+                raise ValueError('GDAL >= 2.1 required to delete nodata values.')
+            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)

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/raster/base.py

@@ -0,0 +1,75 @@
+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,
+                )

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/raster/source.py

@@ -0,0 +1,467 @@
+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_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_text
+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):
+            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_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:
+            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.name.startswith(VSI_FILESYSTEM_BASE_PATH)
+
+    @property
+    def name(self):
+        """
+        Return the name of this raster. Corresponds to filename
+        for file-based rasters.
+        """
+        return force_text(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 transform(self, srid, driver=None, name=None, resampling='NearestNeighbour',
+                  max_error=0.0):
+        """
+        Return a copy of this raster reprojected into the given SRID.
+        """
+        # Convert the resampling algorithm name into an algorithm id
+        algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling]
+
+        # Instantiate target spatial reference system
+        target_srs = SpatialReference(srid)
+
+        # 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': 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.
+        """
+        if not capi.get_ds_info:
+            raise ValueError('GDAL ≥ 2.1 is required for using the info property.')
+        return capi.get_ds_info(self.ptr, None).decode()

venv/lib/python3.7/site-packages/django/contrib/gis/gdal/srs.py

@@ -0,0 +1,335 @@
+"""
+  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 django.contrib.gis.gdal.base import GDALBase
+from django.contrib.gis.gdal.error import SRSException
+from django.contrib.gis.gdal.prototypes import srs as capi
+from django.utils.encoding import force_bytes, force_text
+
+
+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'):
+        """
+        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.4 string, and/or a projection "well known" shorthand
+        string (one of 'WGS84', 'WGS72', 'NAD27', 'NAD83').
+        """
+
+        if srs_type == 'wkt':
+            self.ptr = capi.new_srs(c_char_p(b''))
+            self.import_wkt(srs_input)
+            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
+
+        # 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
+        >>> print(srs['UNIT|AUTHORITY']) # For the units authority, have to use the pipe symbole.
+        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))
+
+    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_text(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.4 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.4 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)