From c078021555dcda6b12363b5b0646eac5332a0d86 Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Fri, 19 Jun 2015 13:56:50 +0100 Subject: [PATCH] Refs #24840 -- Added GDALRaster Warp and transform methods Thanks to Tim Graham for the review. --- django/contrib/gis/gdal/prototypes/raster.py | 3 + django/contrib/gis/gdal/raster/const.py | 11 ++ django/contrib/gis/gdal/raster/source.py | 116 ++++++++++++++++++- docs/ref/contrib/gis/gdal.txt | 86 +++++++++++++- docs/releases/1.9.txt | 9 ++ docs/spelling_wordlist | 1 + tests/gis_tests/gdal_tests/test_raster.py | 116 ++++++++++++++++++- 7 files changed, 337 insertions(+), 5 deletions(-) diff --git a/django/contrib/gis/gdal/prototypes/raster.py b/django/contrib/gis/gdal/prototypes/raster.py index 15b8f64235..12167b4774 100644 --- a/django/contrib/gis/gdal/prototypes/raster.py +++ b/django/contrib/gis/gdal/prototypes/raster.py @@ -69,3 +69,6 @@ get_band_maximum = double_output(std_call('GDALGetRasterMaximum'), [c_void_p, PO reproject_image = void_output(std_call('GDALReprojectImage'), [c_void_p, c_char_p, c_void_p, c_char_p, c_int, c_double, c_double, c_void_p, c_void_p, c_void_p] ) +auto_create_warped_vrt = voidptr_output(std_call('GDALAutoCreateWarpedVRT'), + [c_void_p, c_char_p, c_char_p, c_int, c_double, c_void_p] +) diff --git a/django/contrib/gis/gdal/raster/const.py b/django/contrib/gis/gdal/raster/const.py index 719d1bf0b2..dfe43673b5 100644 --- a/django/contrib/gis/gdal/raster/const.py +++ b/django/contrib/gis/gdal/raster/const.py @@ -32,3 +32,14 @@ GDAL_TO_CTYPES = [ None, c_byte, c_uint16, c_int16, c_uint32, c_int32, c_float, c_double, None, None, None, None ] + +# List of resampling algorithms that can be used to warp a GDALRaster. +GDAL_RESAMPLE_ALGORITHMS = { + 'NearestNeighbour': 0, + 'Bilinear': 1, + 'Cubic': 2, + 'CubicSpline': 3, + 'Lanczos': 4, + 'Average': 5, + 'Mode': 6, +} diff --git a/django/contrib/gis/gdal/raster/source.py b/django/contrib/gis/gdal/raster/source.py index 34f941f40d..e56b442fee 100644 --- a/django/contrib/gis/gdal/raster/source.py +++ b/django/contrib/gis/gdal/raster/source.py @@ -1,12 +1,13 @@ import json import os -from ctypes import addressof, byref, c_double +from ctypes import addressof, byref, c_double, c_void_p 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.raster.const import GDAL_RESAMPLE_ALGORITHMS from django.contrib.gis.gdal.srs import SpatialReference, SRSException from django.contrib.gis.geometry.regex import json_regex from django.utils import six @@ -123,6 +124,9 @@ class GDALRaster(GDALBase): if 'skew' in ds_input: self.skew.x, self.skew.y = ds_input['skew'] + elif isinstance(ds_input, c_void_p): + # Instantiate the object using an existing pointer to a gdal raster. + self._ptr = ds_input else: raise GDALException('Invalid data source input type: "{}".'.format(type(ds_input))) @@ -278,3 +282,113 @@ class GDALRaster(GDALBase): for idx in range(1, capi.get_ds_raster_count(self._ptr) + 1): bands.append(GDALBand(self, idx)) return bands + + def warp(self, ds_input, resampling='NearestNeighbour', max_error=0.0): + """ + Returns a warped GDALRaster with the given input characteristics. + + The input is expected to be a dictionary containing the parameters + of the target raster. Allowed values are width, height, SRID, origin, + scale, skew, datatype, driver, and name (filename). + + By default, the warp functions keeps all parameters equal to the values + of the original source raster. For the name of the target raster, the + name of the source raster will be used and appended with + _copy. + source_driver_name. + + In addition, the resampling algorithm can be specified with the "resampling" + input parameter. The default is NearestNeighbor. For a list of all options + consult the GDAL_RESAMPLE_ALGORITHMS constant. + """ + # Get the parameters defining the geotransform, srid, and size of the raster + if 'width' not in ds_input: + ds_input['width'] = self.width + + if 'height' not in ds_input: + ds_input['height'] = self.height + + if 'srid' not in ds_input: + ds_input['srid'] = self.srs.srid + + if 'origin' not in ds_input: + ds_input['origin'] = self.origin + + if 'scale' not in ds_input: + ds_input['scale'] = self.scale + + if 'skew' not in ds_input: + ds_input['skew'] = self.skew + + # Get the driver, name, and datatype of the target raster + if 'driver' not in ds_input: + ds_input['driver'] = self.driver.name + + if 'name' not in ds_input: + ds_input['name'] = self.name + '_copy.' + self.driver.name + + if 'datatype' not in ds_input: + ds_input['datatype'] = self.bands[0].datatype() + + # Set the number of bands + ds_input['nr_of_bands'] = len(self.bands) + + # Create target raster + target = GDALRaster(ds_input, write=True) + + # Copy nodata values to warped raster + for index, band in enumerate(self.bands): + target.bands[index].nodata_value = band.nodata_value + + # Select resampling algorithm + algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling] + + # Reproject image + capi.reproject_image( + self._ptr, self.srs.wkt.encode(), + target._ptr, target.srs.wkt.encode(), + algorithm, 0.0, max_error, + c_void_p(), c_void_p(), c_void_p() + ) + + # Make sure all data is written to file + target._flush() + + return target + + def transform(self, srid, driver=None, name=None, resampling='NearestNeighbour', + max_error=0.0): + """ + Returns a copy of this raster reprojected into the given SRID. + """ + # Convert the resampling algorithm name into an algorithm id + algorithm = GDAL_RESAMPLE_ALGORITHMS[resampling] + + # Instantiate target spatial reference system + target_srs = SpatialReference(srid) + + # Create warped virtual dataset in the target reference system + target = capi.auto_create_warped_vrt( + self._ptr, self.srs.wkt.encode(), target_srs.wkt.encode(), + algorithm, max_error, c_void_p() + ) + target = GDALRaster(target) + + # Construct the target warp dictionary from the virtual raster + data = { + 'srid': srid, + 'width': target.width, + 'height': target.height, + 'origin': [target.origin.x, target.origin.y], + 'scale': [target.scale.x, target.scale.y], + 'skew': [target.skew.x, target.skew.y], + } + + # Set the driver and filepath if provided + if driver: + data['driver'] = driver + + if name: + data['name'] = name + + # Warp the raster into new srid + return self.warp(data, resampling=resampling, max_error=max_error) diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt index 25f3ecb153..4b68054a75 100644 --- a/docs/ref/contrib/gis/gdal.txt +++ b/docs/ref/contrib/gis/gdal.txt @@ -1119,7 +1119,7 @@ blue. 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. If the input is raw data, the parameters ``width``, - ``heigth``, and ``srid`` are required. The following example shows how rasters + ``height``, and ``srid`` are required. The following example shows how rasters can be created from different input sources (using the sample data from the GeoDjango tests, see also the :ref:`gdal_sample_data` section):: @@ -1288,6 +1288,89 @@ blue. >>> rst.bands[1].data() array([[ 2., 3.]], dtype=float32) + .. method:: warp(ds_input, resampling='NearestNeighbour', max_error=0.0) + + .. versionadded:: 1.9 + + Returns a warped version of this raster. + + The warping parameters can be specified through the ``ds_input`` + argument. The use of ``ds_input`` is analogous to the corresponding + argument of the class constructor. It is a dictionary with the + characteristics of the target raster. Allowed dictionary key values are + width, height, SRID, origin, scale, skew, datatype, driver, and name + (filename). + + By default, the warp functions keeps most parameters equal to the + values of the original source raster, so only parameters that should be + changed need to be specified. Note that this includes the driver, so + for file-based rasters the warp function will create a new raster on + disk. + + The only parameter that is set differently from the source raster is the + name. The default value of the the raster name is the name of the source + raster appended with ``'_copy' + source_driver_name``. For file-based + rasters it is recommended to provide the file path of the target raster. + + The resampling algorithm used for warping can be specified with the + ``resampling`` argument. The default is ``NearestNeighbor``, and the + other allowed values are ``Bilinear``, ``Cubic``, ``CubicSpline``, + ``Lanczos``, ``Average``, and ``Mode``. + + The ``max_error`` argument can be used to specify the maximum error + measured in input pixels that is allowed in approximating the + transformation. The default is 0.0 for exact calculations. + + For users familiar with ``GDAL``, this function has a similar + functionality to the ``gdalwarp`` command-line utility. + + For example, the warp function can be used for aggregating a raster to + the double of its original pixel scale: + + >>> rst = GDALRaster({ + ... "width": 6, "height": 6, "srid": 3086, + ... "origin": [500000, 400000], + ... "scale": [100, -100], + ... "bands": [{"data": range(36), "nodata_value": 99}] + ... }) + >>> target = rst.warp({"scale": [200, -200], "width": 3, "height": 3}) + >>> target.bands[0].data() + array([[ 7., 9., 11.], + [ 19., 21., 23.], + [ 31., 33., 35.]], dtype=float32) + + .. method:: transform(srid, driver=None, name=None, resampling='NearestNeighbour', max_error=0.0) + + .. versionadded:: 1.9 + + Returns a transformed version of this raster with the specified SRID. + + This function transforms the current raster into a new spatial reference + system that can be specified with an ``srid``. It calculates the bounds + and scale of the current raster in the new spatial reference system and + warps the raster using the :attr:`~GDALRaster.warp` function. + + By default, the driver of the source raster is used and the name of the + raster is the original name appended with + ``'_copy' + source_driver_name``. A different driver or name can be + specified with the ``driver`` and ``name`` arguments. + + The default resampling algorithm is ``NearestNeighbour`` but can be + changed using the ``resampling`` argument. The default maximum allowed + error for resampling is 0.0 and can be changed using the ``max_error`` + argument. Consult the :attr:`~GDALRaster.warp` documentation for detail + on those arguments. + + >>> rst = GDALRaster({ + ... "width": 6, "height": 6, "srid": 3086, + ... "origin": [500000, 400000], + ... "scale": [100, -100], + ... "bands": [{"data": range(36), "nodata_value": 99}] + ... }) + >>> target = rst.transform(4326) + >>> target.origin + [-82.98492744885776, 27.601924753080144] + ``GDALBand`` ------------ @@ -1385,7 +1468,6 @@ blue. [ 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 8f3c3bce4e..77d7c09ed3 100644 --- a/docs/releases/1.9.txt +++ b/docs/releases/1.9.txt @@ -177,6 +177,15 @@ Minor features It supports automatic spatial index creation and reprojection when saving a model. It does not yet support spatial querying. +* The new :meth:`GDALRaster.warp() ` + method allows warping a raster by specifying target raster properties such as + origin, width, height, or pixel size (amongst others). + +* The new :meth:`GDALRaster.transform() + ` method allows transforming a + raster into a different spatial reference system by specifying a target + ``srid``. + :mod:`django.contrib.messages` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/spelling_wordlist b/docs/spelling_wordlist index da91c0c25e..293aa6f816 100644 --- a/docs/spelling_wordlist +++ b/docs/spelling_wordlist @@ -660,6 +660,7 @@ repo reportable reprojection reraise +resampling reST reStructuredText reupload diff --git a/tests/gis_tests/gdal_tests/test_raster.py b/tests/gis_tests/gdal_tests/test_raster.py index 1f4847959b..ed1104ed1d 100644 --- a/tests/gis_tests/gdal_tests/test_raster.py +++ b/tests/gis_tests/gdal_tests/test_raster.py @@ -139,8 +139,8 @@ class GDALRasterTests(unittest.TestCase): 'skew': (self.rs.skew.x, self.rs.skew.y), 'bands': [{ 'data': self.rs.bands[0].data(), - 'nodata_value': self.rs.bands[0].nodata_value - }] + 'nodata_value': self.rs.bands[0].nodata_value, + }], }) # Reload newly created raster from file @@ -155,6 +155,118 @@ class GDALRasterTests(unittest.TestCase): else: self.assertEqual(restored_raster.bands[0].data(), self.rs.bands[0].data()) + def test_raster_warp(self): + # Create in memory raster + source = GDALRaster({ + 'datatype': 1, + 'driver': 'MEM', + 'name': 'sourceraster', + 'width': 4, + 'height': 4, + 'nr_of_bands': 1, + 'srid': 3086, + 'origin': (500000, 400000), + 'scale': (100, -100), + 'skew': (0, 0), + 'bands': [{ + 'data': range(16), + 'nodata_value': 255, + }], + }) + + # Test altering the scale, width, and height of a raster + data = { + 'scale': [200, -200], + 'width': 2, + 'height': 2, + } + target = source.warp(data) + self.assertEqual(target.width, data['width']) + self.assertEqual(target.height, data['height']) + self.assertEqual(target.scale, data['scale']) + self.assertEqual(target.bands[0].datatype(), source.bands[0].datatype()) + self.assertEqual(target.name, 'sourceraster_copy.MEM') + result = target.bands[0].data() + if numpy: + result = result.flatten().tolist() + self.assertEqual(result, [5, 7, 13, 15]) + + # Test altering the name and datatype (to float) + data = { + 'name': '/path/to/targetraster.tif', + 'datatype': 6, + } + target = source.warp(data) + self.assertEqual(target.bands[0].datatype(), 6) + self.assertEqual(target.name, '/path/to/targetraster.tif') + self.assertEqual(target.driver.name, 'MEM') + result = target.bands[0].data() + if numpy: + result = result.flatten().tolist() + self.assertEqual( + result, + [0.0, 1.0, 2.0, 3.0, + 4.0, 5.0, 6.0, 7.0, + 8.0, 9.0, 10.0, 11.0, + 12.0, 13.0, 14.0, 15.0] + ) + + def test_raster_transform(self): + # Prepare tempfile and nodata value + rstfile = tempfile.NamedTemporaryFile(suffix='.tif') + ndv = 99 + + # Create in file based raster + source = GDALRaster({ + 'datatype': 1, + 'driver': 'tif', + 'name': rstfile.name, + 'width': 5, + 'height': 5, + 'nr_of_bands': 1, + 'srid': 4326, + 'origin': (-5, 5), + 'scale': (2, -2), + 'skew': (0, 0), + 'bands': [{ + 'data': range(25), + 'nodata_value': ndv, + }], + }) + + # Transform raster into srid 4326. + target = source.transform(3086) + + # Reload data from disk + target = GDALRaster(target.name) + + self.assertEqual(target.srs.srid, 3086) + self.assertEqual(target.width, 7) + self.assertEqual(target.height, 7) + self.assertEqual(target.bands[0].datatype(), source.bands[0].datatype()) + self.assertEqual(target.origin, [9124842.791079799, 1589911.6476407414]) + self.assertEqual(target.scale, [223824.82664250192, -223824.82664250192]) + self.assertEqual(target.skew, [0, 0]) + + result = target.bands[0].data() + if numpy: + result = result.flatten().tolist() + + # The reprojection of a raster that spans over a large area + # skews the data matrix and might introduce nodata values. + self.assertEqual( + result, + [ + ndv, ndv, ndv, ndv, 4, ndv, ndv, + ndv, ndv, 2, 3, 9, ndv, ndv, + ndv, 1, 2, 8, 13, 19, ndv, + 0, 6, 6, 12, 18, 18, 24, + ndv, 10, 11, 16, 22, 23, ndv, + ndv, ndv, 15, 21, 22, ndv, ndv, + ndv, ndv, 20, ndv, ndv, ndv, ndv, + ] + ) + @unittest.skipUnless(HAS_GDAL, "GDAL is required") class GDALBandTests(unittest.TestCase):