Refs #24840 -- Added GDALRaster Warp and transform methods

Thanks to Tim Graham for the review.

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]
+)

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,
+}

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)

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
 ========
 

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() <django.contrib.gis.gdal.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()
+  <django.contrib.gis.gdal.GDALRaster.transform>` method allows transforming a
+  raster into a different spatial reference system by specifying a target
+  ``srid``.
+
 :mod:`django.contrib.messages`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

docs/spelling_wordlist

@@ -660,6 +660,7 @@ repo
 reportable
 reprojection
 reraise
+resampling
 reST
 reStructuredText
 reupload

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):