Added RasterSource/GDALBand GDAL objects

Based on Daniel Wiesmann's raster branch. Thanks Daniel Wiesmann
and Tim Graham for the reviews. Refs #23804.

MANIFEST.in

@@ -26,6 +26,7 @@ recursive-include django/contrib/formtools/tests/templates *
 recursive-include django/contrib/formtools/tests/wizard/wizardtests/templates *
 recursive-include django/contrib/flatpages/fixtures *
 recursive-include django/contrib/flatpages/tests/templates *
+recursive-include django/contrib/gis/gdal/tests/data *
 recursive-include django/contrib/gis/static *
 recursive-include django/contrib/gis/templates *
 recursive-include django/contrib/gis/tests/data *

django/contrib/gis/gdal/__init__.py

@@ -47,6 +47,7 @@ try:
     from django.contrib.gis.gdal.driver import Driver  # NOQA
     from django.contrib.gis.gdal.datasource import DataSource  # NOQA
     from django.contrib.gis.gdal.libgdal import gdal_version, gdal_full_version, GDAL_VERSION  # NOQA
+    from django.contrib.gis.gdal.raster.source import GDALRaster  # NOQA
     from django.contrib.gis.gdal.srs import SpatialReference, CoordTransform  # NOQA
     from django.contrib.gis.gdal.geometries import OGRGeometry  # NOQA
     HAS_GDAL = True

django/contrib/gis/gdal/prototypes/raster.py

@@ -54,8 +54,8 @@ get_band_ds = voidptr_output(lgdal.GDALGetBandDataset, [c_void_p])
 get_band_datatype = int_output(lgdal.GDALGetRasterDataType, [c_void_p])
 get_band_nodata_value = double_output(lgdal.GDALGetRasterNoDataValue, [c_void_p, POINTER(c_int)])
 set_band_nodata_value = void_output(lgdal.GDALSetRasterNoDataValue, [c_void_p, c_double])
-get_band_minimum = double_output(lgdal.GDALGetRasterMinimum, [c_void_p])
-get_band_maximum = double_output(lgdal.GDALGetRasterMaximum, [c_void_p])
+get_band_minimum = double_output(lgdal.GDALGetRasterMinimum, [c_void_p, POINTER(c_int)])
+get_band_maximum = double_output(lgdal.GDALGetRasterMaximum, [c_void_p, POINTER(c_int)])
 
 ### Reprojection routine ###
 reproject_image = void_output(lgdal.GDALReprojectImage, [c_void_p, c_char_p, c_void_p, c_char_p,

django/contrib/gis/gdal/raster/band.py

+from ctypes import byref, c_int
+
+from django.contrib.gis.gdal.base import GDALBase
+from django.contrib.gis.gdal.prototypes import raster as capi
+from django.utils.encoding import force_text
+
+from .const import GDAL_PIXEL_TYPES
+
+
+class GDALBand(GDALBase):
+    """
+    Wraps 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)
+
+    @property
+    def description(self):
+        """
+        Returns 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)
+
+    def datatype(self, as_string=False):
+        """
+        Returns the GDAL Pixel Datatype for this band.
+        """
+        dtype = capi.get_band_datatype(self.ptr)
+        if as_string:
+            dtype = GDAL_PIXEL_TYPES[dtype]
+        return dtype
+
+    @property
+    def min(self):
+        """
+        Returns the minimum pixel value for this band.
+        """
+        return capi.get_band_minimum(self.ptr, byref(c_int()))
+
+    @property
+    def max(self):
+        """
+        Returns the maximum pixel value for this band.
+        """
+        return capi.get_band_maximum(self.ptr, byref(c_int()))
+
+    @property
+    def nodata_value(self):
+        """
+        Returns the nodata value for this band, or None if it isn't set.
+        """
+        nodata_exists = c_int()
+        value = capi.get_band_nodata_value(self.ptr, nodata_exists)
+        return value if nodata_exists else None

django/contrib/gis/gdal/raster/const.py

+"""
+GDAL - Constant definitions
+"""
+
+# See http://www.gdal.org/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4
+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
+}

django/contrib/gis/gdal/raster/source.py

+from ctypes import addressof, byref, c_double
+import os
+
+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.prototypes import raster as capi
+from django.contrib.gis.gdal.raster.band import GDALBand
+from django.contrib.gis.gdal.srs import SpatialReference, SRSException
+from django.utils import six
+from django.utils.six.moves import range
+from django.utils.encoding import (force_bytes, force_text,
+    python_2_unicode_compatible)
+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]]
+        list.__init__(self, [x, y])
+        self._raster = raster
+        self._prop = prop
+
+    @property
+    def x(self):
+        return self[0]
+
+    @property
+    def y(self):
+        return self[1]
+
+
+@python_2_unicode_compatible
+class GDALRaster(GDALBase):
+    """
+    Wraps a raster GDAL Data Source object.
+    """
+    def __init__(self, ds_input, write=False):
+        self._write = 1 if write else 0
+        Driver.ensure_registered()
+
+        # If input is a valid file path, try setting file as source.
+        if isinstance(ds_input, six.string_types):
+            if os.path.exists(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))
+            else:
+                raise GDALException('Unable to read raster source input "{}"'.format(ds_input))
+        else:
+            raise GDALException('Invalid data source input type: "{}".'.format(type(ds_input)))
+
+    def __del__(self):
+        if self._ptr and capi:
+            capi.close_ds(self._ptr)
+
+    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))
+
+    @property
+    def name(self):
+        return force_text(capi.get_ds_description(self.ptr))
+
+    @cached_property
+    def driver(self):
+        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):
+        """
+        Returns the Spatial Reference used in this GDALRaster.
+        """
+        try:
+            wkt = capi.get_ds_projection_ref(self.ptr)
+            return SpatialReference(wkt, srs_type='wkt')
+        except SRSException:
+            return None
+
+    @cached_property
+    def geotransform(self):
+        """
+        Returns the geotransform of the data source.
+        Returns 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 tuple(gtf)
+
+    @property
+    def origin(self):
+        return TransformPoint(self, 'origin')
+
+    @property
+    def scale(self):
+        return TransformPoint(self, 'scale')
+
+    @property
+    def skew(self):
+        return TransformPoint(self, 'skew')
+
+    @property
+    def extent(self):
+        """
+        Returns 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
+
+    @cached_property
+    def bands(self):
+        bands = []
+        for idx in range(1, capi.get_ds_raster_count(self.ptr) + 1):
+            bands.append(GDALBand(self, idx))
+        return bands

django/contrib/gis/gdal/srs.py

@@ -34,7 +34,7 @@ from django.contrib.gis.gdal.error import SRSException
 from django.contrib.gis.gdal.prototypes import srs as capi
 
 from django.utils import six
-from django.utils.encoding import force_bytes
+from django.utils.encoding import force_bytes, force_text
 
 
 #### Spatial Reference class. ####
@@ -46,16 +46,19 @@ class SpatialReference(GDALBase):
     """
 
     #### Python 'magic' routines ####
-    def __init__(self, srs_input=''):
+    def __init__(self, srs_input='', srs_type='user'):
         """
         Creates 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').
         """
-        srs_type = 'user'
 
-        if isinstance(srs_input, six.string_types):
+        if srs_type == 'wkt':
+            self.ptr = capi.new_srs(c_char_p(b''))
+            self.import_wkt(srs_input)
+            return
+        elif isinstance(srs_input, six.string_types):
             # Encoding to ASCII if unicode passed in.
             if isinstance(srs_input, six.text_type):
                 srs_input = srs_input.encode('ascii')
@@ -232,7 +235,7 @@ class SpatialReference(GDALBase):
         elif self.geographic:
             units, name = capi.angular_units(self.ptr, byref(c_char_p()))
         if name is not None:
-            name.decode()
+            name = force_text(name)
         return (units, name)
 
     #### Spheroid/Ellipsoid Properties ####

django/contrib/gis/gdal/tests/test_raster.py

+"""
+gdalinfo django/contrib/gis/gdal/tests/data/raster.tif:
+
+Driver: GTiff/GeoTIFF
+Files: django/contrib/gis/gdal/tests/data/raster.tif
+Size is 163, 174
+Coordinate System is:
+PROJCS["NAD83 / Florida GDL Albers",
+    GEOGCS["NAD83",
+        DATUM["North_American_Datum_1983",
+            SPHEROID["GRS 1980",6378137,298.2572221010002,
+                AUTHORITY["EPSG","7019"]],
+            TOWGS84[0,0,0,0,0,0,0],
+            AUTHORITY["EPSG","6269"]],
+        PRIMEM["Greenwich",0],
+        UNIT["degree",0.0174532925199433],
+        AUTHORITY["EPSG","4269"]],
+    PROJECTION["Albers_Conic_Equal_Area"],
+    PARAMETER["standard_parallel_1",24],
+    PARAMETER["standard_parallel_2",31.5],
+    PARAMETER["latitude_of_center",24],
+    PARAMETER["longitude_of_center",-84],
+    PARAMETER["false_easting",400000],
+    PARAMETER["false_northing",0],
+    UNIT["metre",1,
+        AUTHORITY["EPSG","9001"]],
+    AUTHORITY["EPSG","3086"]]
+Origin = (511700.468070655711927,435103.377123198588379)
+Pixel Size = (100.000000000000000,-100.000000000000000)
+Metadata:
+  AREA_OR_POINT=Area
+Image Structure Metadata:
+  INTERLEAVE=BAND
+Corner Coordinates:
+Upper Left  (  511700.468,  435103.377) ( 82d51'46.16"W, 27d55' 1.53"N)
+Lower Left  (  511700.468,  417703.377) ( 82d51'52.04"W, 27d45'37.50"N)
+Upper Right (  528000.468,  435103.377) ( 82d41'48.81"W, 27d54'56.30"N)
+Lower Right (  528000.468,  417703.377) ( 82d41'55.54"W, 27d45'32.28"N)
+Center      (  519850.468,  426403.377) ( 82d46'50.64"W, 27d50'16.99"N)
+Band 1 Block=163x50 Type=Byte, ColorInterp=Gray
+  NoData Value=15
+"""
+import os
+import unittest
+
+from django.contrib.gis.gdal import HAS_GDAL
+from django.utils import six
+from django.utils._os import upath
+
+if HAS_GDAL:
+    from django.contrib.gis.gdal import GDALRaster
+    from django.contrib.gis.gdal.raster.band import GDALBand
+
+
+@unittest.skipUnless(HAS_GDAL, "GDAL is required")
+class GDALRasterTests(unittest.TestCase):
+    """
+    Test a GDALRaster instance created from a file (GeoTiff).
+    """
+    def setUp(self):
+        self.rs_path = os.path.join(os.path.dirname(upath(__file__)),
+                                    'data/raster.tif')
+        self.rs = GDALRaster(self.rs_path)
+
+    def test_rs_name_repr(self):
+        self.assertEqual(self.rs_path, self.rs.name)
+        six.assertRegex(self, repr(self.rs), "<Raster object at 0x\w+>")
+
+    def test_rs_driver(self):
+        self.assertEqual(self.rs.driver.name, 'GTiff')
+
+    def test_rs_size(self):
+        self.assertEqual(self.rs.width, 163)
+        self.assertEqual(self.rs.height, 174)
+
+    def test_rs_srs(self):
+        self.assertEqual(self.rs.srs.srid, 3086)
+        self.assertEqual(self.rs.srs.units, (1.0, 'metre'))
+
+    def test_geotransform_and_friends(self):
+        self.assertEqual(self.rs.geotransform,
+            (511700.4680706557, 100.0, 0.0, 435103.3771231986, 0.0, -100.0))
+        self.assertEqual(self.rs.origin, [511700.4680706557, 435103.3771231986])
+        self.assertEqual(self.rs.origin.x, 511700.4680706557)
+        self.assertEqual(self.rs.origin.y, 435103.3771231986)
+        self.assertEqual(self.rs.scale, [100.0, -100.0])
+        self.assertEqual(self.rs.scale.x, 100.0)
+        self.assertEqual(self.rs.scale.y, -100.0)
+        self.assertEqual(self.rs.skew, [0, 0])
+        self.assertEqual(self.rs.skew.x, 0)
+        self.assertEqual(self.rs.skew.y, 0)
+
+    def test_rs_extent(self):
+        self.assertEqual(self.rs.extent,
+            (511700.4680706557, 417703.3771231986, 528000.4680706557, 435103.3771231986))
+
+    def test_rs_bands(self):
+        self.assertEqual(len(self.rs.bands), 1)
+        self.assertIsInstance(self.rs.bands[0], GDALBand)
+
+
+@unittest.skipUnless(HAS_GDAL, "GDAL is required")
+class GDALBandTests(unittest.TestCase):
+    def setUp(self):
+        rs_path = os.path.join(os.path.dirname(upath(__file__)),
+                               'data/raster.tif')
+        rs = GDALRaster(rs_path)
+        self.band = rs.bands[0]
+
+    def test_band_data(self):
+        self.assertEqual(self.band.width, 163)
+        self.assertEqual(self.band.height, 174)
+        self.assertEqual(self.band.description, '')
+        self.assertEqual(self.band.datatype(), 1)
+        self.assertEqual(self.band.datatype(as_string=True), 'GDT_Byte')
+        self.assertEqual(self.band.min, 0)
+        self.assertEqual(self.band.max, 255)
+        self.assertEqual(self.band.nodata_value, 15)

docs/ref/contrib/gis/gdal.txt

@@ -18,8 +18,8 @@ of vector spatial data.
 .. note::
 
      Although the module is named ``gdal``, GeoDjango only supports
-     some of the capabilities of OGR.  Thus, none of GDAL's features
-     with respect to raster (image) data are supported at this time.
+     some of the capabilities of OGR.  Thus, GDAL's features with respect to
+     raster (image) data are minimally supported (read-only) at this time.
 
 __ http://www.gdal.org/
 __ http://www.gdal.org/ogr/
@@ -1081,6 +1081,140 @@ the same coordinate transformation repeatedly on different geometries::
     ...     geom = feat.geom # getting clone of feature geometry
     ...     geom.transform(ct) # transforming
 
+.. _raster-data-source-objects:
+
+Raster Data Objects
+===================
+
+.. versionadded:: 1.8
+
+``GDALRaster``
+----------------
+
+:class:`GDALRaster` is a wrapper for the GDAL raster source object that
+supports reading data from a variety of GDAL-supported geospatial file
+formats and data sources using a simple, consistent interface.  Each
+data source is represented by a :class:`GDALRaster` object which contains
+one or more layers of data named bands.  Each band, represented by a
+:class:`GDALBand` object, contains georeferenced image data. For exemple, an RGB
+image is represented as three bands: one for red, one for green, and one for
+blue.
+
+.. class:: GDALRaster(ds_input)
+
+    The constructor for ``GDALRaster`` accepts a single parameter: the path of
+    the file you want to read.
+
+    .. attribute:: name
+
+        The name of the source which is equivalent to the input file path.
+
+    .. attribute:: driver
+
+        The name of the GDAL driver used to handle the input file. For example,
+        ``GTiff`` for a ``GeoTiff`` file. See also the `GDAL Raster Formats`__
+        list.
+
+        __ http://www.gdal.org/formats_list.html
+
+    .. attribute:: width
+
+       The width of the source in pixels (X-axis).
+
+    .. attribute:: height
+
+        The height of the source in pixels (Y-axis).
+
+    .. attribute:: srs
+
+        The spatial reference system of the source, as a
+        :class:`SpatialReference` instance.
+
+    .. attribute:: geotransform
+
+        The affine transformation matrix used to georeference the source, as a
+        tuple of six coefficients which map pixel/line coordinates into
+        georeferenced space using the following relationship::
+
+            Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
+            Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
+
+        The same values can be retrieved by accessing the :attr:`origin`
+        (indices 0 and 3), :attr:`scale` (indices 1 and 5) and :attr:`skew`
+        (indices 2 and 4) properties.
+
+    .. attribute:: origin
+
+        Coordinates of the top left origin of the raster in the spatial
+        reference system of the source, as a point object with ``x`` and ``y``
+        members.
+
+    .. attribute:: scale
+
+        Pixel width and height used for georeferencing the raster, as a as a
+        point object with ``x`` and ``y``  members. See :attr:`geotransform`
+        for more information.
+
+    .. attribute:: skew
+
+        Skew coefficients used to georeference the raster, as a point object
+        with ``x`` and ``y``  members. In case of north up images, these
+        coefficients are both ``0``.
+
+    .. attribute:: extent
+
+        Extent (boundary values) of the raster source, as a 4-tuple
+        ``(xmin, ymin, xmax, ymax)`` in the spatial reference system of the
+        source.
+
+    .. attribute:: bands
+
+        List of all bands of the source, as :class:`GDALBand` instances.
+
+``GDALBand``
+------------
+
+.. class:: GDALBand
+
+    ``GDALBand`` instances are not created explicitely, but rather obtained
+    from a :class:`GDALRaster` object, through its :attr:`~GDALRaster.bands`
+    attribute.
+
+    .. attribute:: description
+
+        The name or description of the band, if any.
+
+    .. attribute:: width
+
+        The width of the band in pixels (X-axis).
+
+    .. attribute:: height
+
+        The height of the band in pixels (Y-axis).
+
+    .. attribute:: min
+
+        The minimum pixel value of the band (excluding the "no data" value).
+
+    .. attribute:: max
+
+        The maximum pixel value of the band (excluding the "no data" value).
+
+    .. attribute:: nodata_value
+
+        The "no data" value for a band is generally a special marker value used
+        to mark pixels that are not valid data. Such pixels should generally not
+        be displayed, nor contribute to analysis operations.
+
+    .. method:: datatype([as_string=False])
+
+        The data type contained in the band, as an integer constant between 0
+        (Unknown) and 11. If ``as_string`` is ``True``, the data type is
+        returned as a string with the following possible values:
+        ``GDT_Unknown``, ``GDT_Byte``, ``GDT_UInt16``, ``GDT_Int16``,
+        ``GDT_UInt32``, ``GDT_Int32``, ``GDT_Float32``, ``GDT_Float64``,
+        ``GDT_CInt16``, ``GDT_CInt32``, ``GDT_CFloat32``, and ``GDT_CFloat64``.
+
 Settings
 ========
 

docs/releases/1.8.txt

@@ -166,6 +166,9 @@ Minor features
   ``SELECT InitSpatialMetaData`` initialization commands are now automatically
   run by :djadmin:`migrate`.
 
+* The GDAL interface now supports retrieving properties of
+  :ref:`raster (image) data file <raster-data-source-objects>`.
+
 * Compatibility shims for ``SpatialRefSys`` and ``GeometryColumns`` changed in
   Django 1.2 have been removed.