diff --git a/MANIFEST.in b/MANIFEST.in index d96430ce67..47f90d28f8 100644 --- a/MANIFEST.in +++ b/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 * diff --git a/django/contrib/gis/gdal/__init__.py b/django/contrib/gis/gdal/__init__.py index 1258dfcdae..ff6ed7efcf 100644 --- a/django/contrib/gis/gdal/__init__.py +++ b/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 diff --git a/django/contrib/gis/gdal/prototypes/raster.py b/django/contrib/gis/gdal/prototypes/raster.py index 180952276f..db32ee1e0b 100644 --- a/django/contrib/gis/gdal/prototypes/raster.py +++ b/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, diff --git a/django/contrib/gis/gdal/raster/__init__.py b/django/contrib/gis/gdal/raster/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/django/contrib/gis/gdal/raster/band.py b/django/contrib/gis/gdal/raster/band.py new file mode 100644 index 0000000000..2c7aaab533 --- /dev/null +++ b/django/contrib/gis/gdal/raster/band.py @@ -0,0 +1,69 @@ +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 diff --git a/django/contrib/gis/gdal/raster/const.py b/django/contrib/gis/gdal/raster/const.py new file mode 100644 index 0000000000..c98ac771cb --- /dev/null +++ b/django/contrib/gis/gdal/raster/const.py @@ -0,0 +1,19 @@ +""" +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 +} diff --git a/django/contrib/gis/gdal/raster/source.py b/django/contrib/gis/gdal/raster/source.py new file mode 100644 index 0000000000..e89e70d2e8 --- /dev/null +++ b/django/contrib/gis/gdal/raster/source.py @@ -0,0 +1,155 @@ +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 '' % 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 diff --git a/django/contrib/gis/gdal/srs.py b/django/contrib/gis/gdal/srs.py index 548a3cd0ea..8749d079b3 100644 --- a/django/contrib/gis/gdal/srs.py +++ b/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 #### diff --git a/django/contrib/gis/gdal/tests/data/raster.tif b/django/contrib/gis/gdal/tests/data/raster.tif new file mode 100644 index 0000000000..8a50ae24e6 Binary files /dev/null and b/django/contrib/gis/gdal/tests/data/raster.tif differ diff --git a/django/contrib/gis/gdal/tests/test_raster.py b/django/contrib/gis/gdal/tests/test_raster.py new file mode 100644 index 0000000000..4dd6042ee5 --- /dev/null +++ b/django/contrib/gis/gdal/tests/test_raster.py @@ -0,0 +1,118 @@ +""" +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), "") + + 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) diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt index ce3cdc79d5..6491e7cb76 100644 --- a/docs/ref/contrib/gis/gdal.txt +++ b/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 ======== diff --git a/docs/releases/1.8.txt b/docs/releases/1.8.txt index 77609ed7b5..3dc075a7c8 100644 --- a/docs/releases/1.8.txt +++ b/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 `. + * Compatibility shims for ``SpatialRefSys`` and ``GeometryColumns`` changed in Django 1.2 have been removed.