From f269c1d6f6dcc22c0a781f3223c6da0a4483b06e Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Fri, 13 Mar 2015 18:49:02 +0000 Subject: [PATCH] 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 | 1 + django/contrib/gis/gdal/raster/band.py | 103 +++++++-- django/contrib/gis/gdal/raster/const.py | 12 ++ django/contrib/gis/gdal/raster/source.py | 160 ++++++++++++-- docs/ref/contrib/gis/gdal.txt | 197 ++++++++++++++++-- docs/releases/1.9.txt | 5 +- tests/gis_tests/data/__init__.py | 0 tests/gis_tests/data/rasters/__init__.py | 0 .../data => data/rasters}/raster.tif | Bin 29158 -> 35486 bytes tests/gis_tests/data/rasters/textrasters.py | 12 ++ tests/gis_tests/gdal_tests/test_raster.py | 167 ++++++++++++++- 11 files changed, 600 insertions(+), 57 deletions(-) create mode 100644 tests/gis_tests/data/__init__.py create mode 100644 tests/gis_tests/data/rasters/__init__.py rename tests/gis_tests/{gdal_tests/data => data/rasters}/raster.tif (81%) create mode 100644 tests/gis_tests/data/rasters/textrasters.py diff --git a/django/contrib/gis/gdal/prototypes/raster.py b/django/contrib/gis/gdal/prototypes/raster.py index 08ec9df961..a79e2a9db3 100644 --- a/django/contrib/gis/gdal/prototypes/raster.py +++ b/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]) diff --git a/django/contrib/gis/gdal/raster/band.py b/django/contrib/gis/gdal/raster/band.py index 2c7aaab533..3f6c0a4140 100644 --- a/django/contrib/gis/gdal/raster/band.py +++ b/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) - if as_string: - dtype = GDAL_PIXEL_TYPES[dtype] - return dtype + return self.width * self.height @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() diff --git a/django/contrib/gis/gdal/raster/const.py b/django/contrib/gis/gdal/raster/const.py index c98ac771cb..bdf858afaf 100644 --- a/django/contrib/gis/gdal/raster/const.py +++ b/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 +] diff --git a/django/contrib/gis/gdal/raster/source.py b/django/contrib/gis/gdal/raster/source.py index bb72db49ee..06f6b5b03f 100644 --- a/django/contrib/gis/gdal/raster/source.py +++ b/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): - 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: + if not os.path.exists(ds_input): 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 '' % hex(addressof(self.ptr)) + return '' % 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)) - return tuple(gtf) + capi.get_ds_geotransform(self._ptr, byref(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 diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt index 365b6a99c9..ef2b83b48c 100644 --- a/docs/ref/contrib/gis/gdal.txt +++ b/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 - raster (image) data are minimally supported (read-only) at this time. + some of the capabilities of OGR and GDAL's raster features 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 - the file you want to read. + For raster data there is no difference between a raster instance and its + 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, - ``GTiff`` for a ``GeoTiff`` file. See also the `GDAL Raster Formats`__ - list. + The name of the GDAL driver used to handle the input file. For ``GDALRaster``\s created + from a file, the driver type is detected automatically. The creation of rasters from + 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 - :class:`SpatialReference` instance. + The spatial reference system of the raster, as a + :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 ======== diff --git a/docs/releases/1.9.txt b/docs/releases/1.9.txt index 4375ba2a0d..1f9caa1145 100644 --- a/docs/releases/1.9.txt +++ b/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 ` from raw data. + Setters for raster properties such as projection or pixel values have + been added. :mod:`django.contrib.messages` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/tests/gis_tests/data/__init__.py b/tests/gis_tests/data/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/gis_tests/data/rasters/__init__.py b/tests/gis_tests/data/rasters/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/gis_tests/gdal_tests/data/raster.tif b/tests/gis_tests/data/rasters/raster.tif similarity index 81% rename from tests/gis_tests/gdal_tests/data/raster.tif rename to tests/gis_tests/data/rasters/raster.tif index 8a50ae24e604bc56e9df86912a410c23960c86b3..ddc7605703246def0c3fd01390dcb5e763cceb38 100644 GIT binary patch delta 1277 zcmajbu}cDB9LDjtOQfK|z&}7l8lu4>nj#tsDjXad9UL1PPBT06Of&ONkwA{XL7EyG zgbErQ8XFuMC8Q}z!YLXW8hhS1e*B;gZu?Fj9-jBNv5RxfzqvUpp6y%zOKNPoIk~R& z(#?2S>!zJEJ-jkV3sJ623BM5H6Mb!siwHWjC&Z04C4A@|^sS|dFxo~RUut3;Js)*G z^cm+Owb6W9a<)?2zkVYleKRBdg?5JhJ2d%2H2HJ;?8o+$XuA2urY4&nz^C-urIAk-5CeuKZ z*+P5#+|EnCkyo2LdG;UBwD}87{<6t_!gSr7EJ!n6P@DAvn>L!vE}Ben#AXQXu{lzd zez2%EFBaL4p~*MVVRBm{e8; delta 24 gcmbO?mFd}IMov#pEe53m1_r*#`K7KK{bm#b0A?fzDF6Tf diff --git a/tests/gis_tests/data/rasters/textrasters.py b/tests/gis_tests/data/rasters/textrasters.py new file mode 100644 index 0000000000..b12ce62afc --- /dev/null +++ b/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]}] +} +""" diff --git a/tests/gis_tests/gdal_tests/test_raster.py b/tests/gis_tests/gdal_tests/test_raster.py index 4dd6042ee5..36b4aa7547 100644 --- a/tests/gis_tests/gdal_tests/test_raster.py +++ b/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__)), - 'data/raster.tif') - rs = GDALRaster(rs_path) + self.rs_path = os.path.join(os.path.dirname(upath(__file__)), + '../data/rasters/raster.tif') + 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))