Added write support for GDALRaster

- Instantiation of GDALRaster instances from dict or json data.
- Retrieve and write pixel values in GDALBand objects.
- Support for the GDALFlushCache in gdal C prototypes
- Added private flush method to GDALRaster to make sure all
  data is written to files when file-based rasters are changed.
- Replaced ``ptr`` with ``_ptr`` for internal ptr variable

Refs #23804. Thanks Claude Paroz and Tim Graham for the reviews.

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

@@ -31,6 +31,7 @@ get_driver_description = const_string_output(lgdal.GDALGetDescription, [c_void_p
 create_ds = voidptr_output(lgdal.GDALCreate, [c_void_p, c_char_p, c_int, c_int, c_int, c_int])
 open_ds = voidptr_output(lgdal.GDALOpen, [c_char_p, c_int])
 close_ds = void_output(lgdal.GDALClose, [c_void_p])
+flush_ds = int_output(lgdal.GDALFlushCache, [c_void_p])
 copy_ds = voidptr_output(lgdal.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(lgdal.GDALAddBand, [c_void_p, c_int])

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

@@ -2,9 +2,11 @@ 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.contrib.gis.shortcuts import numpy
+from django.utils import six
 from django.utils.encoding import force_text
 
-from .const import GDAL_PIXEL_TYPES
+from .const import GDAL_PIXEL_TYPES, GDAL_TO_CTYPES
 
 
 class GDALBand(GDALBase):
@@ -13,51 +15,49 @@ class GDALBand(GDALBase):
     """
     def __init__(self, source, index):
         self.source = source
-        self.ptr = capi.get_ds_raster_band(source.ptr, index)
+        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))
+        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)
+        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)
+        return capi.get_band_ysize(self._ptr)
 
-    def datatype(self, as_string=False):
+    @property
+    def pixel_count(self):
         """
-        Returns the GDAL Pixel Datatype for this band.
+        Returns the total number of pixels in this band.
         """
-        dtype = capi.get_band_datatype(self.ptr)
+        return self.width * self.height
-        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()))
+        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()))
+        return capi.get_band_maximum(self._ptr, byref(c_int()))
 
     @property
     def nodata_value(self):
@@ -65,5 +65,80 @@ class GDALBand(GDALBase):
         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)
+        value = capi.get_band_nodata_value(self._ptr, nodata_exists)
         return value if nodata_exists else None
+
+    @nodata_value.setter
+    def nodata_value(self, value):
+        """
+        Sets the nodata value for this band.
+        """
+        if not isinstance(value, (int, float)):
+            raise ValueError('Nodata value must be numeric.')
+        capi.set_band_nodata_value(self._ptr, value)
+        self.source._flush()
+
+    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
+
+    def data(self, data=None, offset=None, size=None, as_memoryview=False):
+        """
+        Reads 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.
+        """
+        if not offset:
+            offset = (0, 0)
+
+        if not size:
+            size = (self.width - offset[0], self.height - offset[1])
+
+        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()] * (size[0] * size[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, six.memoryview, 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), size[0],
+                     size[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:
+                return numpy.frombuffer(
+                    data_array, dtype=numpy.dtype(data_array)).reshape(size)
+            else:
+                return list(data_array)
+        else:
+            self.source._flush()

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

@@ -1,6 +1,9 @@
 """
 GDAL - Constant definitions
 """
+from ctypes import (
+    c_byte, c_double, c_float, c_int16, c_int32, c_uint16, c_uint32,
+)
 
 # See http://www.gdal.org/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4
 GDAL_PIXEL_TYPES = {
@@ -17,3 +20,12 @@ GDAL_PIXEL_TYPES = {
     10: 'GDT_CFloat32',  # Complex Float32
     11: 'GDT_CFloat64',  # Complex Float64
 }
+
+# Lookup values to convert GDAL pixel type indices into ctypes objects.
+# The GDAL band-io works with ctypes arrays to hold data to be written
+# or to hold the space for data to be read into. The lookup below helps
+# selecting the right ctypes object for a given gdal pixel type.
+GDAL_TO_CTYPES = [
+    None, c_byte, c_uint16, c_int16, c_uint32, c_int32,
+    c_float, c_double, None, None, None, None
+]

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

@@ -1,3 +1,4 @@
+import json
 import os
 from ctypes import addressof, byref, c_double
 
@@ -7,6 +8,7 @@ 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.contrib.gis.geometry.regex import json_regex
 from django.utils import six
 from django.utils.encoding import (
     force_bytes, force_text, python_2_unicode_compatible,
@@ -33,10 +35,22 @@ class TransformPoint(list):
     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
+
 
 @python_2_unicode_compatible
 class GDALRaster(GDALBase):
@@ -47,17 +61,64 @@ class GDALRaster(GDALBase):
         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, six.string_types) 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, six.string_types):
-            if os.path.exists(ds_input):
+            if not 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))
+            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, 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.')
+
+            # 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),
+                None
+            )
+
+            # Set band data if provided
+            for i, band_input in enumerate(ds_input.get('bands', [])):
+                self.bands[i].data(band_input['data'])
+                if 'nodata_value' in band_input:
+                    self.bands[i].nodata_value = band_input['nodata_value']
+
+            # Set SRID, default to 0 (this assures SRS is always instanciated)
+            self.srs = ds_input.get('srid', 0)
+
+            # 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']
         else:
             raise GDALException('Invalid data source input type: "{}".'.format(type(ds_input)))
 
@@ -72,15 +133,34 @@ class GDALRaster(GDALBase):
         """
         Short-hand representation because WKB may be very large.
         """
-        return '<Raster object at %s>' % hex(addressof(self.ptr))
+        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 name(self):
-        return force_text(capi.get_ds_description(self.ptr))
+        """
+        Returns 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):
-        ds_driver = capi.get_ds_driver(self.ptr)
+        """
+        Returns the GDAL Driver used for this raster.
+        """
+        ds_driver = capi.get_ds_driver(self._ptr)
         return Driver(ds_driver)
 
     @property
@@ -88,48 +168,85 @@ class GDALRaster(GDALBase):
         """
         Width (X axis) in pixels.
         """
-        return capi.get_ds_xsize(self.ptr)
+        return capi.get_ds_xsize(self._ptr)
 
     @property
     def height(self):
         """
         Height (Y axis) in pixels.
         """
-        return capi.get_ds_ysize(self.ptr)
+        return capi.get_ds_ysize(self._ptr)
 
     @property
     def srs(self):
         """
-        Returns the Spatial Reference used in this GDALRaster.
+        Returns the SpatialReference used in this GDALRaster.
         """
         try:
-            wkt = capi.get_ds_projection_ref(self.ptr)
+            wkt = capi.get_ds_projection_ref(self._ptr)
+            if not wkt:
+                return None
             return SpatialReference(wkt, srs_type='wkt')
         except SRSException:
             return None
 
-    @cached_property
+    @srs.setter
+    def srs(self, value):
+        """
+        Sets 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, six.integer_types + six.string_types):
+            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 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).
+        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))
+        capi.get_ds_geotransform(self._ptr, byref(gtf))
-        return tuple(gtf)
+        return list(gtf)
+
+    @geotransform.setter
+    def geotransform(self, values):
+        "Sets the geotransform for the data source."
+        if sum([isinstance(x, (int, float)) for x in values]) != 6:
+            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
@@ -150,7 +267,10 @@ class GDALRaster(GDALBase):
 
     @cached_property
     def bands(self):
+        """
+        Returns the bands of this raster as a list of GDALBand instances.
+        """
         bands = []
-        for idx in range(1, capi.get_ds_raster_count(self.ptr) + 1):
+        for idx in range(1, capi.get_ds_raster_count(self._ptr) + 1):
             bands.append(GDALBand(self, idx))
         return bands

docs/ref/contrib/gis/gdal.txt

@@ -13,13 +13,13 @@ formats.
 
 GeoDjango provides a high-level Python interface for some of the
 capabilities of OGR, including the reading and coordinate transformation
-of vector spatial data.
+of vector spatial data and minimal support for GDAL's features with respect
+to raster (image) data.
 
 .. note::
 
      Although the module is named ``gdal``, GeoDjango only supports
-     some of the capabilities of OGR.  Thus, GDAL's features with respect to
+     some of the capabilities of OGR and GDAL's raster features at this time.
-     raster (image) data are minimally supported (read-only) at this time.
 
 __ http://www.gdal.org/
 __ http://www.gdal.org/ogr/
@@ -27,6 +27,8 @@ __ http://www.gdal.org/ogr/
 Overview
 ========
 
+.. _gdal_sample_data:
+
 Sample Data
 -----------
 
@@ -37,6 +39,7 @@ have any data of your own to use, GeoDjango tests contain a number of
 simple data sets that you can use for testing. You can download them here::
 
     $ wget https://raw.githubusercontent.com/django/django/master/tests/gis_tests/data/cities/cities.{shp,prj,shx,dbf}
+    $ wget https://raw.githubusercontent.com/django/django/master/tests/gis_tests/data/rasters/raster.tif
 
 Vector Data Source Objects
 ==========================
@@ -1101,35 +1104,106 @@ one or more layers of data named bands.  Each band, represented by a
 image is represented as three bands: one for red, one for green, and one for
 blue.
 
-.. class:: GDALRaster(ds_input)
+.. note::
 
-    The constructor for ``GDALRaster`` accepts a single parameter: the path of
+    For raster data there is no difference between a raster instance and its
-    the file you want to read.
+    data source. Unlike for the Geometry objects, :class:`GDALRaster` objects are
+    always a data source. Temporary rasters can be instantiated in memory
+    using the corresponding driver, but they will be of the same class as file-based
+    raster sources.
+
+.. class:: GDALRaster(ds_input, write=False)
+
+    The constructor for ``GDALRaster`` accepts two parameters. The first parameter
+    defines the raster source, it is either a path to a file or spatial data with
+    values defining the properties of a new raster (such as size and name). If the
+    input is a file path, the second parameter specifies if the raster should
+    be opened with write access. The following example shows how rasters can be
+    created from different input sources (using the sample data from the GeoDjango
+    tests, see the :ref:`gdal_sample_data` section)::
+
+        >>> from django.contrib.gis.gdal.raster.source import GDALRaster
+        >>> rst = GDALRaster('/path/to/your/raster.tif', write=False)
+        >>> rst.name
+        '/path/to/your/raster.tif'
+        >>> rst.width, rst.height            # This file has 163 x 174 pixels
+        (163, 174)
+        >>> rst = GDALRaster({'srid': 4326, 'width': 1, 'height': 2, 'datatype': 1
+        ...                   'bands': [{'data': [0, 1]}]}) # Creates in-memory raster
+        >>> rst.srs.srid
+        4326
+        >>> rst.width, rst.height
+        (1, 2)
+        >>> rst.bands[0].data()
+        array([[0, 1]], dtype=int8)
+
+    .. versionchanged:: 1.9
+
+        ``GDALRaster`` objects can now be instantiated directly from raw data.
+        Setters have been added for the following properties: ``srs``,
+        ``geotransform``, ``origin``, ``scale``, and ``skew``.
 
     .. attribute:: name
 
-        The name of the source which is equivalent to the input file path.
+        The name of the source which is equivalent to the input file path or the name
+        provided upon instantiation.
+
+            >>> GDALRaster({'width': 10, 'height': 10, 'name': 'myraster'}).name
+            'myraster'
 
     .. attribute:: driver
 
-        The name of the GDAL driver used to handle the input file. For example,
+        The name of the GDAL driver used to handle the input file. For ``GDALRaster``\s created
-        ``GTiff`` for a ``GeoTiff`` file. See also the `GDAL Raster Formats`__
+        from a file, the driver type is detected automatically. The creation of rasters from
-        list.
+        scratch is a in-memory raster by default (``'MEM'``), but can be altered as
+        needed. For instance, use ``GTiff`` for a ``GeoTiff`` file. For a list of file types,
+        see also the `GDAL Raster Formats`__ list.
 
         __ http://www.gdal.org/formats_list.html
 
+        An in-memory raster is created through the following example:
+
+            >>> GDALRaster({'width': 10, 'height': 10}).driver.name
+            'MEM'
+
+        A file based GeoTiff raster is created through the following example:
+
+            >>> import tempfile
+            >>> rstfile = tempfile.NamedTemporaryFile(suffix='.tif')
+            >>> rst = GDALRaster({'driver': 'GTiff', 'name': rstfile.name,
+            ...                   'width': 255, 'height': 255, 'nr_of_bands': 1})
+            >>> rst.name
+            '/tmp/tmp7x9H4J.tif'           # The exact filename will be different on your computer
+            >>> rst.driver.name
+            'GTiff'
+
     .. attribute:: width
 
-       The width of the source in pixels (X-axis).
+        The width of the source in pixels (X-axis).
+
+            >>> GDALRaster({'width': 10, 'height': 20}).width
+            10
 
     .. attribute:: height
 
         The height of the source in pixels (Y-axis).
 
+            >>> GDALRaster({'width': 10, 'height': 20}).height
+            20
+
     .. attribute:: srs
 
-        The spatial reference system of the source, as a
+        The spatial reference system of the raster, as a
-        :class:`SpatialReference` instance.
+        :class:`SpatialReference` instance. The SRS can be changed by
+        setting it to an other :class:`SpatialReference` or providing any input
+        that is accepted by the :class:`SpatialReference` constructor.
+
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.srs
+            None
+            >>> rst.srs = 4326
+            >>> rst.srs.srid
+            4326
 
     .. attribute:: geotransform
 
@@ -1144,34 +1218,75 @@ blue.
         (indices 0 and 3), :attr:`scale` (indices 1 and 5) and :attr:`skew`
         (indices 2 and 4) properties.
 
+        The default is ``[0.0, 1.0, 0.0, 0.0, 0.0, -1.0]``.
+
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.geotransform
+            [0.0, 1.0, 0.0, 0.0, 0.0, -1.0]
+
     .. 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.
 
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.origin
+            [0.0, 0.0]
+            >>> rst.origin.x = 1
+            >>> rst.origin
+            [1.0, 0.0]
+
     .. 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.
 
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.scale
+            [1.0, -1.0]
+            >>> rst.scale.x = 2
+            >>> rst.scale
+            [2.0, -1.0]
+
     .. 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``.
 
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.skew
+            [0.0, 0.0]
+            >>> rst.skew.x = 3
+            >>> rst.skew
+            [3.0, 0.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.
 
+            >>> rst = GDALRaster({'width': 10, 'height': 20})
+            >>> rst.extent
+            (0.0, -20.0, 10.0, 0.0)
+            >>> rst.origin.x = 100
+            >>> rst.extent
+            (100.0, -20.0, 110.0, 0.0)
+
     .. attribute:: bands
 
         List of all bands of the source, as :class:`GDALBand` instances.
 
+            >>> rst = GDALRaster({"width": 1, "height": 2, "bands": [{"data": [0, 1]},
+            ...                                                      {"data": [2, 3]}]})
+            >>> len(rst.bands)
+            2
+            >>> rst.bands[1].data()
+            array([[ 2.,  3.]], dtype=float32)
+
 ``GDALBand``
 ------------
 
@@ -1179,7 +1294,7 @@ blue.
 
     ``GDALBand`` instances are not created explicitly, but rather obtained
     from a :class:`GDALRaster` object, through its :attr:`~GDALRaster.bands`
-    attribute.
+    attribute. The GDALBands contain the actual pixel values of the raster.
 
     .. attribute:: description
 
@@ -1193,6 +1308,12 @@ blue.
 
         The height of the band in pixels (Y-axis).
 
+    .. attribute:: pixel_count
+
+        .. versionadded:: 1.9
+
+        The total number of pixels in this band. Is equal to ``width * height``.
+
     .. attribute:: min
 
         The minimum pixel value of the band (excluding the "no data" value).
@@ -1207,6 +1328,10 @@ blue.
         to mark pixels that are not valid data. Such pixels should generally not
         be displayed, nor contribute to analysis operations.
 
+        .. versionchanged:: 1.9
+
+            This property can now be set as well.
+
     .. method:: datatype([as_string=False])
 
         The data type contained in the band, as an integer constant between 0
@@ -1216,6 +1341,50 @@ blue.
         ``GDT_UInt32``, ``GDT_Int32``, ``GDT_Float32``, ``GDT_Float64``,
         ``GDT_CInt16``, ``GDT_CInt32``, ``GDT_CFloat32``, and ``GDT_CFloat64``.
 
+    .. method:: data(data=None, offset=None, size=None)
+
+        .. versionadded:: 1.9
+
+        The accessor to the pixel values of the ``GDALBand``. Returns the complete
+        data array if no parameters are provided. A subset of the pixel array can
+        be requested by specifying an offset and block size as tuples.
+
+        If NumPy is available, the data is returned as NumPy array. For performance
+        reasons, it is highly recommended to use NumPy.
+
+        Data is written to the ``GDALBand`` if the ``data`` parameter is provided.
+        The input can be of one of the following types - packed string, buffer, list,
+        array, and NumPy array. The number of items in the input must correspond to the
+        total number of pixels in the band, or to the number of pixels for a specific
+        block of pixel values if the ``offset`` and ``size`` parameters are provided.
+
+        For example:
+
+            >>> rst = GDALRaster({'width': 4, 'height': 4, 'datatype': 1, 'nr_of_bands': 1})
+            >>> bnd = rst.bands[0]
+            >>> bnd.data(range(16))
+            >>> bnd.data()
+            array([[ 0,  1,  2,  3],
+                   [ 4,  5,  6,  7],
+                   [ 8,  9, 10, 11],
+                   [12, 13, 14, 15]], dtype=int8)
+            >>> bnd.data(offset=(1,1), size=(2,2))
+            array([[ 5,  6],
+                   [ 9, 10]], dtype=int8)
+            >>> bnd.data(data=[-1, -2, -3, -4], offset=(1,1), size=(2,2))
+            >>> bnd.data()
+            array([[ 0,  1,  2,  3],
+                   [ 4, -1, -2,  7],
+                   [ 8, -3, -4, 11],
+                   [12, 13, 14, 15]], dtype=int8)
+            >>> bnd.data(data='\x9d\xa8\xb3\xbe', offset=(1,1), size=(2,2))
+            >>> bnd.data()
+            array([[  0,   1,   2,   3],
+                   [  4, -99, -88,   7],
+                   [  8, -77, -66,  11],
+                   [ 12,  13,  14,  15]], dtype=int8)
+
+
 Settings
 ========
 

docs/releases/1.9.txt

@@ -54,7 +54,10 @@ Minor features
 :mod:`django.contrib.gis`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-* ...
+* The GDAL interface now supports instantiating file-based and in-memory
+  :ref:`GDALRaster objects <raster-data-source-objects>` from raw data.
+  Setters for raster properties such as projection or pixel values have
+  been added.
 
 :mod:`django.contrib.messages`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

tests/gis_tests/data/rasters/textrasters.py

@@ -0,0 +1,12 @@
+JSON_RASTER = """{
+    "srid": 4326,
+    "origin": [0, 0],
+    "scale": [1, 1],
+    "skew": [0, 0],
+    "width": 5,
+    "height": 5,
+    "nr_of_bands": 1,
+    "bands": [{"data": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
+                        15, 16, 17, 18, 19, 20, 21, 22, 23, 24]}]
+}
+"""

tests/gis_tests/gdal_tests/test_raster.py

@@ -1,8 +1,8 @@
 """
-gdalinfo django/contrib/gis/gdal/tests/data/raster.tif:
+gdalinfo tests/gis_tests/data/rasters/raster.tif:
 
 Driver: GTiff/GeoTIFF
-Files: django/contrib/gis/gdal/tests/data/raster.tif
+Files: tests/gis_tests/data/rasters/raster.tif
 Size is 163, 174
 Coordinate System is:
 PROJCS["NAD83 / Florida GDL Albers",
@@ -41,12 +41,18 @@ Band 1 Block=163x50 Type=Byte, ColorInterp=Gray
   NoData Value=15
 """
 import os
+import struct
+import tempfile
 import unittest
 
 from django.contrib.gis.gdal import HAS_GDAL
+from django.contrib.gis.gdal.error import GDALException
+from django.contrib.gis.shortcuts import numpy
 from django.utils import six
 from django.utils._os import upath
 
+from ..data.rasters.textrasters import JSON_RASTER
+
 if HAS_GDAL:
     from django.contrib.gis.gdal import GDALRaster
     from django.contrib.gis.gdal.raster.band import GDALBand
@@ -59,7 +65,7 @@ class GDALRasterTests(unittest.TestCase):
     """
     def setUp(self):
         self.rs_path = os.path.join(os.path.dirname(upath(__file__)),
-                                    'data/raster.tif')
+                                    '../data/rasters/raster.tif')
         self.rs = GDALRaster(self.rs_path)
 
     def test_rs_name_repr(self):
@@ -78,8 +84,9 @@ class GDALRasterTests(unittest.TestCase):
         self.assertEqual(self.rs.srs.units, (1.0, 'metre'))
 
     def test_geotransform_and_friends(self):
+        # Assert correct values for file based raster
         self.assertEqual(self.rs.geotransform,
-            (511700.4680706557, 100.0, 0.0, 435103.3771231986, 0.0, -100.0))
+            [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)
@@ -89,22 +96,72 @@ class GDALRasterTests(unittest.TestCase):
         self.assertEqual(self.rs.skew, [0, 0])
         self.assertEqual(self.rs.skew.x, 0)
         self.assertEqual(self.rs.skew.y, 0)
+        # Create in-memory rasters and change gtvalues
+        rsmem = GDALRaster(JSON_RASTER)
+        rsmem.geotransform = range(6)
+        self.assertEqual(rsmem.geotransform, [float(x) for x in range(6)])
+        self.assertEqual(rsmem.origin, [0, 3])
+        self.assertEqual(rsmem.origin.x, 0)
+        self.assertEqual(rsmem.origin.y, 3)
+        self.assertEqual(rsmem.scale, [1, 5])
+        self.assertEqual(rsmem.scale.x, 1)
+        self.assertEqual(rsmem.scale.y, 5)
+        self.assertEqual(rsmem.skew, [2, 4])
+        self.assertEqual(rsmem.skew.x, 2)
+        self.assertEqual(rsmem.skew.y, 4)
+        self.assertEqual(rsmem.width, 5)
+        self.assertEqual(rsmem.height, 5)
 
     def test_rs_extent(self):
         self.assertEqual(self.rs.extent,
-            (511700.4680706557, 417703.3771231986, 528000.4680706557, 435103.3771231986))
+            (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)
 
+    def test_file_based_raster_creation(self):
+        # Prepare tempfile
+        rstfile = tempfile.NamedTemporaryFile(suffix='.tif')
+
+        # Create file-based raster from scratch
+        GDALRaster({
+            'datatype': self.rs.bands[0].datatype(),
+            'driver': 'tif',
+            'name': rstfile.name,
+            'width': 163,
+            'height': 174,
+            'nr_of_bands': 1,
+            'srid': self.rs.srs.wkt,
+            'origin': (self.rs.origin.x, self.rs.origin.y),
+            'scale': (self.rs.scale.x, self.rs.scale.y),
+            'skew': (self.rs.skew.x, self.rs.skew.y),
+            'bands': [{
+                'data': self.rs.bands[0].data(),
+                'nodata_value': self.rs.bands[0].nodata_value
+            }]
+        })
+
+        # Reload newly created raster from file
+        restored_raster = GDALRaster(rstfile.name)
+        self.assertEqual(restored_raster.srs.wkt, self.rs.srs.wkt)
+        self.assertEqual(restored_raster.geotransform, self.rs.geotransform)
+        if numpy:
+            numpy.testing.assert_equal(
+                restored_raster.bands[0].data(),
+                self.rs.bands[0].data()
+            )
+        else:
+            self.assertEqual(restored_raster.bands[0].data(), self.rs.bands[0].data())
+
 
 @unittest.skipUnless(HAS_GDAL, "GDAL is required")
 class GDALBandTests(unittest.TestCase):
     def setUp(self):
-        rs_path = os.path.join(os.path.dirname(upath(__file__)),
+        self.rs_path = os.path.join(os.path.dirname(upath(__file__)),
-                               'data/raster.tif')
+                               '../data/rasters/raster.tif')
-        rs = GDALRaster(rs_path)
+        rs = GDALRaster(self.rs_path)
         self.band = rs.bands[0]
 
     def test_band_data(self):
@@ -116,3 +173,97 @@ class GDALBandTests(unittest.TestCase):
         self.assertEqual(self.band.min, 0)
         self.assertEqual(self.band.max, 255)
         self.assertEqual(self.band.nodata_value, 15)
+
+    def test_read_mode_error(self):
+        # Open raster in read mode
+        rs = GDALRaster(self.rs_path, write=False)
+        band = rs.bands[0]
+
+        # Setting attributes in write mode raises exception in the _flush method
+        self.assertRaises(GDALException, setattr, band, 'nodata_value', 10)
+
+    def test_band_data_setters(self):
+        # Create in-memory raster and get band
+        rsmem = GDALRaster({
+            'datatype': 1,
+            'driver': 'MEM',
+            'name': 'mem_rst',
+            'width': 10,
+            'height': 10,
+            'nr_of_bands': 1
+        })
+        bandmem = rsmem.bands[0]
+
+        # Set nodata value
+        bandmem.nodata_value = 99
+        self.assertEqual(bandmem.nodata_value, 99)
+
+        # Set data for entire dataset
+        bandmem.data(range(100))
+        if numpy:
+            numpy.testing.assert_equal(bandmem.data(), numpy.arange(100).reshape(10, 10))
+        else:
+            self.assertEqual(bandmem.data(), range(100))
+
+        # Prepare data for setting values in subsequent tests
+        block = range(100, 104)
+        packed_block = struct.pack('<' + 'B B B B', *block)
+
+        # Set data from list
+        bandmem.data(block, (1, 1), (2, 2))
+        result = bandmem.data(offset=(1, 1), size=(2, 2))
+        if numpy:
+            numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))
+        else:
+            self.assertEqual(result, block)
+
+        # Set data from packed block
+        bandmem.data(packed_block, (1, 1), (2, 2))
+        result = bandmem.data(offset=(1, 1), size=(2, 2))
+        if numpy:
+            numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))
+        else:
+            self.assertEqual(result, block)
+
+        # Set data from bytes
+        bandmem.data(bytes(packed_block), (1, 1), (2, 2))
+        result = bandmem.data(offset=(1, 1), size=(2, 2))
+        if numpy:
+            numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))
+        else:
+            self.assertEqual(result, block)
+
+        # Set data from bytearray
+        bandmem.data(bytearray(packed_block), (1, 1), (2, 2))
+        result = bandmem.data(offset=(1, 1), size=(2, 2))
+        if numpy:
+            numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))
+        else:
+            self.assertEqual(result, block)
+
+        # Set data from memoryview
+        bandmem.data(six.memoryview(packed_block), (1, 1), (2, 2))
+        result = bandmem.data(offset=(1, 1), size=(2, 2))
+        if numpy:
+            numpy.testing.assert_equal(result, numpy.array(block).reshape(2, 2))
+        else:
+            self.assertEqual(result, block)
+
+        # Set data from numpy array
+        if numpy:
+            bandmem.data(numpy.array(block, dtype='int8').reshape(2, 2), (1, 1), (2, 2))
+            numpy.testing.assert_equal(
+                bandmem.data(offset=(1, 1), size=(2, 2)),
+                numpy.array(block).reshape(2, 2)
+            )
+
+        # Test json input data
+        rsmemjson = GDALRaster(JSON_RASTER)
+        bandmemjson = rsmemjson.bands[0]
+        if numpy:
+            numpy.testing.assert_equal(
+                bandmemjson.data(),
+                numpy.array(range(25)).reshape(5, 5)
+            )
+        else:
+            self.assertEqual(bandmemjson.data(), range(25))