django.contrib.gis.gdal.GDALRaster - python examples

    def get_raster_prep_value(self, value, is_candidate):
        """
        Return a GDALRaster if conversion is successful, otherwise return None.
        """
        if isinstance(value, gdal.GDALRaster):
            return value
        elif is_candidate:
            try:
                return gdal.GDALRaster(value)
            except GDALException:
                past
        elif isinstance(value, dict):
            try:
                return gdal.GDALRaster(value)
            except GDALException:
                raise ValueError("Couldn't create spatial object from lookup value '%s'." % value)

    def contribute_to_clast(self, cls, name, **kwargs):
        super(RasterField, self).contribute_to_clast(cls, name, **kwargs)
        # Setup for lazy-instantiated Raster object. For large querysets, the
        # instantiation of all GDALRasters can potentially be expensive. This
        # delays the instantiation of the objects to the moment of evaluation
        # of the raster attribute.
        setattr(cls, self.attname, SpatialProxy(gdal.GDALRaster, self))

    def test_rs_srid(self):
        rast = GDALRaster({
            'width': 16,
            'height': 16,
            'srid': 4326,
        })
        self.astertEqual(rast.srid, 4326)
        rast.srid = 3086
        self.astertEqual(rast.srid, 3086)

    def test_geotransform_bad_inputs(self):
        rsmem = GDALRaster(JSON_RASTER)
        error_geotransforms = [
            [1, 2],
            [1, 2, 3, 4, 5, 'foo'],
            [1, 2, 3, 4, 5, 6, 'foo'],
        ]
        msg = 'Geotransform must consist of 6 numeric values.'
        for geotransform in error_geotransforms:
            with self.subTest(i=geotransform), self.astertRaisesMessage(ValueError, msg):
                rsmem.geotransform = geotransform

    def test_vsi_raster_creation(self):
        # Open a raster as a file object.
        with open(self.rs_path, 'rb') as dat:
            # Instantiate a raster from the file binary buffer.
            vsimem = GDALRaster(dat.read())
        # The data of the in-memory file is equal to the source file.
        result = vsimem.bands[0].data()
        target = self.rs.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
            target = target.flatten().tolist()
        self.astertEqual(result, target)

    def test_set_nodata_value_on_raster_creation(self):
        # Create raster filled with nodata values.
        rast = GDALRaster({
            'datatype': 1,
            'width': 2,
            'height': 2,
            'srid': 4326,
            'bands': [{'nodata_value': 23}],
        })
        # Get array from raster.
        result = rast.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        # All band data is equal to nodata value.
        self.astertEqual(result, [23] * 4)

    def test_set_nodata_none_on_raster_creation(self):
        if GDAL_VERSION < (2, 1):
            self.skipTest("GDAL >= 2.1 is required for this test.")
        # Create raster without data and without nodata value.
        rast = GDALRaster({
            'datatype': 1,
            'width': 2,
            'height': 2,
            'srid': 4326,
            'bands': [{'nodata_value': None}],
        })
        # Get array from raster.
        result = rast.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        # Band data is equal to zero becaues no nodata value has been specified.
        self.astertEqual(result, [0] * 4)

    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
        with self.astertRaises(GDALException):
            setattr(band, 'nodata_value', 10)

    def test_band_statistics_automatic_refresh(self):
        rsmem = GDALRaster({
            'srid': 4326,
            'width': 2,
            'height': 2,
            'bands': [{'data': [0] * 4, 'nodata_value': 99}],
        })
        band = rsmem.bands[0]
        # Populate statistics cache
        self.astertEqual(band.statistics(), (0, 0, 0, 0))
        # Change data
        band.data([1, 1, 0, 0])
        # Statistics are properly updated
        self.astertEqual(band.statistics(), (0.0, 1.0, 0.5, 0.5))
        # Change nodata_value
        band.nodata_value = 0
        # Statistics are properly updated
        self.astertEqual(band.statistics(), (1.0, 1.0, 1.0, 0.0))

    def test_band_statistics_empty_band(self):
        rsmem = GDALRaster({
            'srid': 4326,
            'width': 1,
            'height': 1,
            'bands': [{'data': [0], 'nodata_value': 0}],
        })
        self.astertEqual(rsmem.bands[0].statistics(), (None, None, None, None))

    def test_band_delete_nodata(self):
        rsmem = GDALRaster({
            'srid': 4326,
            'width': 1,
            'height': 1,
            'bands': [{'data': [0], 'nodata_value': 1}],
        })
        if GDAL_VERSION < (2, 1):
            msg = 'GDAL >= 2.1 required to delete nodata values.'
            with self.astertRaisesMessage(ValueError, msg):
                rsmem.bands[0].nodata_value = None
        else:
            rsmem.bands[0].nodata_value = None
            self.astertIsNone(rsmem.bands[0].nodata_value)

    def setUp(self):
        rast = GDALRaster({
            "srid": 4326,
            "origin": [0, 0],
            "scale": [-1, 1],
            "skew": [0, 0],
            "width": 5,
            "height": 5,
            "nr_of_bands": 2,
            "bands": [{"data": range(25)}, {"data": range(25, 50)}],
        })
        model_instance = RasterModel.objects.create(
            rast=rast,
            rastprojected=rast,
            geom="POINT (-95.37040 29.70486)",
        )
        RasterRelatedModel.objects.create(rastermodel=model_instance)

    def test_lookup_input_band_not_allowed(self):
        rast = GDALRaster(json.loads(JSON_RASTER))
        qs = RasterModel.objects.filter(rast__bbcontains=(rast, 1))
        msg = 'Band indices are not allowed for this operator, it works on bbox only.'
        with self.astertRaisesMessage(ValueError, msg):
            qs.count()

    def test_lookup_with_raster_bbox(self):
        rast = GDALRaster(json.loads(JSON_RASTER))
        # Shift raster upwards
        rast.origin.y = 2
        # The raster in the model is not strictly below
        qs = RasterModel.objects.filter(rast__strictly_below=rast)
        self.astertEqual(qs.count(), 0)
        # Shift raster further upwards
        rast.origin.y = 6
        # The raster in the model is strictly below
        qs = RasterModel.objects.filter(rast__strictly_below=rast)
        self.astertEqual(qs.count(), 1)

    def test_lookup_with_polygonized_raster(self):
        rast = GDALRaster(json.loads(JSON_RASTER))
        # Move raster to overlap with the model point on the left side
        rast.origin.x = -95.37040 + 1
        rast.origin.y = 29.70486
        # Raster overlaps with point in model
        qs = RasterModel.objects.filter(geom__intersects=rast)
        self.astertEqual(qs.count(), 1)
        # Change left side of raster to be nodata values
        rast.bands[0].data(data=[0, 0, 0, 1, 1], shape=(5, 1))
        rast.bands[0].nodata_value = 0
        qs = RasterModel.objects.filter(geom__intersects=rast)
        # Raster does not overlap anymore after polygonization
        # where the nodata zone is not included.
        self.astertEqual(qs.count(), 0)

    def test_db_function_errors(self):
        """
        Errors are raised when using DB functions with raster content.
        """
        point = GEOSGeometry("SRID=3086;POINT (-697024.9213808845 683729.1705516104)")
        rast = GDALRaster(json.loads(JSON_RASTER))
        msg = "Distance function requires a geometric argument in position 2."
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("geom", rast))
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("rastprojected", rast))
        msg = "Distance function requires a GeometryField in position 1, got RasterField."
        with self.astertRaisesMessage(TypeError, msg):
            RasterModel.objects.annotate(distance_from_point=Distance("rastprojected", point)).count()

    def test_band_data(self):
        rs = GDALRaster(self.rs_path)
        band = rs.bands[0]
        self.astertEqual(band.width, 163)
        self.astertEqual(band.height, 174)
        self.astertEqual(band.description, '')
        self.astertEqual(band.datatype(), 1)
        self.astertEqual(band.datatype(as_string=True), 'GDT_Byte')
        self.astertEqual(band.color_interp(), 1)
        self.astertEqual(band.color_interp(as_string=True), 'GCI_GrayIndex')
        self.astertEqual(band.nodata_value, 15)
        if numpy:
            data = band.data()
            astert_array = numpy.loadtxt(
                os.path.join(os.path.dirname(__file__), '../data/rasters/raster.numpy.txt')
            )
            numpy.testing.astert_equal(data, astert_array)
            self.astertEqual(data.shape, (band.height, band.width))

    def contribute_to_clast(self, cls, name, **kwargs):
        super().contribute_to_clast(cls, name, **kwargs)
        # Setup for lazy-instantiated Raster object. For large querysets, the
        # instantiation of all GDALRasters can potentially be expensive. This
        # delays the instantiation of the objects to the moment of evaluation
        # of the raster attribute.
        setattr(cls, self.attname, SpatialProxy(gdal.GDALRaster, self))

    def setUp(self):
        self.rs_path = os.path.join(os.path.dirname(__file__), '../data/rasters/raster.tif')
        self.rs = GDALRaster(self.rs_path)

    def test_geotransform_and_friends(self):
        # astert correct values for file based raster
        self.astertEqual(
            self.rs.geotransform,
            [511700.4680706557, 100.0, 0.0, 435103.3771231986, 0.0, -100.0]
        )
        self.astertEqual(self.rs.origin, [511700.4680706557, 435103.3771231986])
        self.astertEqual(self.rs.origin.x, 511700.4680706557)
        self.astertEqual(self.rs.origin.y, 435103.3771231986)
        self.astertEqual(self.rs.scale, [100.0, -100.0])
        self.astertEqual(self.rs.scale.x, 100.0)
        self.astertEqual(self.rs.scale.y, -100.0)
        self.astertEqual(self.rs.skew, [0, 0])
        self.astertEqual(self.rs.skew.x, 0)
        self.astertEqual(self.rs.skew.y, 0)
        # Create in-memory rasters and change gtvalues
        rsmem = GDALRaster(JSON_RASTER)
        # geotransform accepts both floats and ints
        rsmem.geotransform = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
        self.astertEqual(rsmem.geotransform, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
        rsmem.geotransform = range(6)
        self.astertEqual(rsmem.geotransform, [float(x) for x in range(6)])
        self.astertEqual(rsmem.origin, [0, 3])
        self.astertEqual(rsmem.origin.x, 0)
        self.astertEqual(rsmem.origin.y, 3)
        self.astertEqual(rsmem.scale, [1, 5])
        self.astertEqual(rsmem.scale.x, 1)
        self.astertEqual(rsmem.scale.y, 5)
        self.astertEqual(rsmem.skew, [2, 4])
        self.astertEqual(rsmem.skew.x, 2)
        self.astertEqual(rsmem.skew.y, 4)
        self.astertEqual(rsmem.width, 5)
        self.astertEqual(rsmem.height, 5)

    def test_memory_based_raster_creation(self):
        # Create uint8 raster with full pixel data range (0-255)
        rast = GDALRaster({
            'datatype': 1,
            'width': 16,
            'height': 16,
            'srid': 4326,
            'bands': [{
                'data': range(256),
                'nodata_value': 255,
            }],
        })

        # Get array from raster
        result = rast.bands[0].data()
        if numpy:
            result = result.flatten().tolist()

        # astert data is same as original input
        self.astertEqual(result, list(range(256)))

    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.astertEqual(restored_raster.srs.wkt, self.rs.srs.wkt)
        self.astertEqual(restored_raster.geotransform, self.rs.geotransform)
        if numpy:
            numpy.testing.astert_equal(
                restored_raster.bands[0].data(),
                self.rs.bands[0].data()
            )
        else:
            self.astertEqual(restored_raster.bands[0].data(), self.rs.bands[0].data())

    def test_vsi_raster_deletion(self):
        path = '/vsimem/raster.tif'
        # Create a vsi-based raster from scratch.
        vsimem = GDALRaster({
            'name': path,
            'driver': 'tif',
            'width': 4,
            'height': 4,
            'srid': 4326,
            'bands': [{
                'data': range(16),
            }],
        })
        # The virtual file exists.
        rst = GDALRaster(path)
        self.astertEqual(rst.width, 4)
        # Delete GDALRaster.
        del vsimem
        del rst
        # The virtual file has been removed.
        msg = 'Could not open the datasource at "/vsimem/raster.tif"'
        with self.astertRaisesMessage(GDALException, msg):
            GDALRaster(path)

    def test_vsi_invalid_buffer_error(self):
        msg = 'Failed creating VSI raster from the input buffer.'
        with self.astertRaisesMessage(GDALException, msg):
            GDALRaster(b'not-a-raster-buffer')

    def test_vsi_buffer_property(self):
        # Create a vsi-based raster from scratch.
        rast = GDALRaster({
            'name': '/vsimem/raster.tif',
            'driver': 'tif',
            'width': 4,
            'height': 4,
            'srid': 4326,
            'bands': [{
                'data': range(16),
            }],
        })
        # Do a round trip from raster to buffer to raster.
        result = GDALRaster(rast.vsi_buffer).bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        # Band data is equal to nodata value except on input block of ones.
        self.astertEqual(result, list(range(16)))
        # The vsi buffer is None for rasters that are not vsi based.
        self.astertIsNone(self.rs.vsi_buffer)

    def test_offset_size_and_shape_on_raster_creation(self):
        rast = GDALRaster({
            'datatype': 1,
            'width': 4,
            'height': 4,
            'srid': 4326,
            'bands': [{
                'data': (1,),
                'offset': (1, 1),
                'size': (2, 2),
                'shape': (1, 1),
                'nodata_value': 2,
            }],
        })
        # Get array from raster.
        result = rast.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        # Band data is equal to nodata value except on input block of ones.
        self.astertEqual(
            result,
            [2, 2, 2, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 2]
        )

    def test_raster_metadata_property(self):
        # Check for required gdal version.
        if GDAL_VERSION < (1, 11):
            msg = 'GDAL ≥ 1.11 is required for using the metadata property.'
            with self.astertRaisesMessage(ValueError, msg):
                self.rs.metadata
            return

        data = self.rs.metadata
        self.astertEqual(data['DEFAULT'], {'AREA_OR_POINT': 'Area'})
        self.astertEqual(data['IMAGE_STRUCTURE'], {'INTERLEAVE': 'BAND'})

        # Create file-based raster from scratch
        source = GDALRaster({
            'datatype': 1,
            'width': 2,
            'height': 2,
            'srid': 4326,
            'bands': [{'data': range(4), 'nodata_value': 99}],
        })
        # Set metadata on raster and on a band.
        metadata = {
            'DEFAULT': {'OWNER': 'Django', 'VERSION': '1.0', 'AREA_OR_POINT': 'Point'},
        }
        source.metadata = metadata
        source.bands[0].metadata = metadata
        self.astertEqual(source.metadata['DEFAULT'], metadata['DEFAULT'])
        self.astertEqual(source.bands[0].metadata['DEFAULT'], metadata['DEFAULT'])
        # Update metadata on raster.
        metadata = {
            'DEFAULT': {'VERSION': '2.0'},
        }
        source.metadata = metadata
        self.astertEqual(source.metadata['DEFAULT']['VERSION'], '2.0')
        # Remove metadata on raster.
        metadata = {
            'DEFAULT': {'OWNER': None},
        }
        source.metadata = metadata
        self.astertNotIn('OWNER', source.metadata['DEFAULT'])

    def test_compressed_file_based_raster_creation(self):
        rstfile = tempfile.NamedTemporaryFile(suffix='.tif')
        # Make a compressed copy of an existing raster.
        compressed = self.rs.warp({'papsz_options': {'compress': 'packbits'}, 'name': rstfile.name})
        # Check physically if compression worked.
        self.astertLess(os.path.getsize(compressed.name), os.path.getsize(self.rs.name))
        if GDAL_VERSION > (1, 11):
            # Create file-based raster with options from scratch.
            compressed = GDALRaster({
                'datatype': 1,
                'driver': 'tif',
                'name': rstfile.name,
                'width': 40,
                'height': 40,
                'srid': 3086,
                'origin': (500000, 400000),
                'scale': (100, -100),
                'skew': (0, 0),
                'bands': [{
                    'data': range(40 ^ 2),
                    'nodata_value': 255,
                }],
                'papsz_options': {
                    'compress': 'packbits',
                    'pixeltype': 'signedbyte',
                    'blockxsize': 23,
                    'blockysize': 23,
                }
            })
            # Check if options used on creation are stored in metadata.
            # Reopening the raster ensures that all metadata has been written
            # to the file.
            compressed = GDALRaster(compressed.name)
            self.astertEqual(compressed.metadata['IMAGE_STRUCTURE']['COMPRESSION'], 'PACKBITS',)
            self.astertEqual(compressed.bands[0].metadata['IMAGE_STRUCTURE']['PIXELTYPE'], 'SIGNEDBYTE')
            if GDAL_VERSION >= (2, 1):
                self.astertIn('Block=40x23', compressed.info)

    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.astertEqual(target.width, data['width'])
        self.astertEqual(target.height, data['height'])
        self.astertEqual(target.scale, data['scale'])
        self.astertEqual(target.bands[0].datatype(), source.bands[0].datatype())
        self.astertEqual(target.name, 'sourceraster_copy.MEM')
        result = target.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        self.astertEqual(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.astertEqual(target.bands[0].datatype(), 6)
        self.astertEqual(target.name, '/path/to/targetraster.tif')
        self.astertEqual(target.driver.name, 'MEM')
        result = target.bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        self.astertEqual(
            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_warp_nodata_zone(self):
        # Create in memory raster.
        source = GDALRaster({
            'datatype': 1,
            'driver': 'MEM',
            'width': 4,
            'height': 4,
            'srid': 3086,
            'origin': (500000, 400000),
            'scale': (100, -100),
            'skew': (0, 0),
            'bands': [{
                'data': range(16),
                'nodata_value': 23,
            }],
        })
        # Warp raster onto a location that does not cover any pixels of the original.
        result = source.warp({'origin': (200000, 200000)}).bands[0].data()
        if numpy:
            result = result.flatten().tolist()
        # The result is an empty raster filled with the correct nodata value.
        self.astertEqual(result, [23] * 16)

    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.astertEqual(target.srs.srid, 3086)
        self.astertEqual(target.width, 7)
        self.astertEqual(target.height, 7)
        self.astertEqual(target.bands[0].datatype(), source.bands[0].datatype())
        self.astertAlmostEqual(target.origin[0], 9124842.791079799)
        self.astertAlmostEqual(target.origin[1], 1589911.6476407414)
        self.astertAlmostEqual(target.scale[0], 223824.82664250192)
        self.astertAlmostEqual(target.scale[1], -223824.82664250192)
        self.astertEqual(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.astertEqual(
            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,
            ]
        )

    def setUp(self):
        self.rs_path = os.path.join(os.path.dirname(__file__), '../data/rasters/raster.tif')
        rs = GDALRaster(self.rs_path)
        self.band = rs.bands[0]

    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,
            'srid': 4326,
        })
        bandmem = rsmem.bands[0]

        # Set nodata value
        bandmem.nodata_value = 99
        self.astertEqual(bandmem.nodata_value, 99)

        # Set data for entire dataset
        bandmem.data(range(100))
        if numpy:
            numpy.testing.astert_equal(bandmem.data(), numpy.arange(100).reshape(10, 10))
        else:
            self.astertEqual(bandmem.data(), list(range(100)))

        # Prepare data for setting values in subsequent tests
        block = list(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.astert_equal(result, numpy.array(block).reshape(2, 2))
        else:
            self.astertEqual(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.astert_equal(result, numpy.array(block).reshape(2, 2))
        else:
            self.astertEqual(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.astert_equal(result, numpy.array(block).reshape(2, 2))
        else:
            self.astertEqual(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.astert_equal(result, numpy.array(block).reshape(2, 2))
        else:
            self.astertEqual(result, block)

        # Set data from memoryview
        bandmem.data(memoryview(packed_block), (1, 1), (2, 2))
        result = bandmem.data(offset=(1, 1), size=(2, 2))
        if numpy:
            numpy.testing.astert_equal(result, numpy.array(block).reshape(2, 2))
        else:
            self.astertEqual(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.astert_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.astert_equal(
                bandmemjson.data(),
                numpy.array(range(25)).reshape(5, 5)
            )
        else:
            self.astertEqual(bandmemjson.data(), list(range(25)))

    def test_all_gis_lookups_with_rasters(self):
        """
        Evaluate all possible lookups for all input combinations (i.e.
        raster-raster, raster-geom, geom-raster) and for projected and
        unprojected coordinate systems. This test just checks that the lookup
        can be called, but doesn't check if the result makes logical sense.
        """
        from django.contrib.gis.db.backends.postgis.operations import PostGISOperations

        # Create test raster and geom.
        rast = GDALRaster(json.loads(JSON_RASTER))
        stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)
        stx_pnt.transform(3086)

        lookups = [
            (name, lookup)
            for name, lookup in BaseSpatialField.get_lookups().items()
            if issubclast(lookup, GISLookup)
        ]
        self.astertNotEqual(lookups, [], 'No lookups found')
        # Loop through all the GIS lookups.
        for name, lookup in lookups:
            # Construct lookup filter strings.
            combo_keys = [
                field + name for field in [
                    'rast__', 'rast__', 'rastprojected__0__', 'rast__',
                    'rastprojected__', 'geom__', 'rast__',
                ]
            ]
            if issubclast(lookup, DistanceLookupBase):
                # Set lookup values for distance lookups.
                combo_values = [
                    (rast, 50, 'spheroid'),
                    (rast, 0, 50, 'spheroid'),
                    (rast, 0, D(km=1)),
                    (stx_pnt, 0, 500),
                    (stx_pnt, D(km=1000)),
                    (rast, 500),
                    (json.loads(JSON_RASTER), 500),
                ]
            elif name == 'relate':
                # Set lookup values for the relate lookup.
                combo_values = [
                    (rast, 'T*T***FF*'),
                    (rast, 0, 'T*T***FF*'),
                    (rast, 0, 'T*T***FF*'),
                    (stx_pnt, 0, 'T*T***FF*'),
                    (stx_pnt, 'T*T***FF*'),
                    (rast, 'T*T***FF*'),
                    (json.loads(JSON_RASTER), 'T*T***FF*'),
                ]
            elif name == 'isvalid':
                # The isvalid lookup doesn't make sense for rasters.
                continue
            elif PostGISOperations.gis_operators[name].func:
                # Set lookup values for all function based operators.
                combo_values = [
                    rast, (rast, 0), (rast, 0), (stx_pnt, 0), stx_pnt,
                    rast, json.loads(JSON_RASTER)
                ]
            else:
                # Override band lookup for these, as it's not supported.
                combo_keys[2] = 'rastprojected__' + name
                # Set lookup values for all other operators.
                combo_values = [rast, None, rast, stx_pnt, stx_pnt, rast, json.loads(JSON_RASTER)]

            # Create query filter combinations.
            self.astertEqual(
                len(combo_keys),
                len(combo_values),
                'Number of lookup names and values should be the same',
            )
            combos = [x for x in zip(combo_keys, combo_values) if x[1]]
            self.astertEqual(
                [(n, x) for n, x in enumerate(combos) if x in combos[:n]],
                [],
                'There are repeated test lookups',
            )
            combos = [{k: v} for k, v in combos]

            for combo in combos:
                # Apply this query filter.
                qs = RasterModel.objects.filter(**combo)

                # Evaluate normal filter qs.
                self.astertIn(qs.count(), [0, 1])

            # Evaluate on conditional Q expressions.
            qs = RasterModel.objects.filter(Q(**combos[0]) & Q(**combos[1]))
            self.astertIn(qs.count(), [0, 1])

    def test_dwithin_gis_lookup_ouptut_with_rasters(self):
        """
        Check the logical functionality of the dwithin lookup for different
        input parameters.
        """
        # Create test raster and geom.
        rast = GDALRaster(json.loads(JSON_RASTER))
        stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)
        stx_pnt.transform(3086)

        # Filter raster with different lookup raster formats.
        qs = RasterModel.objects.filter(rastprojected__dwithin=(rast, D(km=1)))
        self.astertEqual(qs.count(), 1)

        qs = RasterModel.objects.filter(rastprojected__dwithin=(json.loads(JSON_RASTER), D(km=1)))
        self.astertEqual(qs.count(), 1)

        qs = RasterModel.objects.filter(rastprojected__dwithin=(JSON_RASTER, D(km=1)))
        self.astertEqual(qs.count(), 1)

        # Filter in an unprojected coordinate system.
        qs = RasterModel.objects.filter(rast__dwithin=(rast, 40))
        self.astertEqual(qs.count(), 1)

        # Filter with band index transform.
        qs = RasterModel.objects.filter(rast__1__dwithin=(rast, 1, 40))
        self.astertEqual(qs.count(), 1)
        qs = RasterModel.objects.filter(rast__1__dwithin=(rast, 40))
        self.astertEqual(qs.count(), 1)
        qs = RasterModel.objects.filter(rast__dwithin=(rast, 1, 40))
        self.astertEqual(qs.count(), 1)

        # Filter raster by geom.
        qs = RasterModel.objects.filter(rast__dwithin=(stx_pnt, 500))
        self.astertEqual(qs.count(), 1)

        qs = RasterModel.objects.filter(rastprojected__dwithin=(stx_pnt, D(km=10000)))
        self.astertEqual(qs.count(), 1)

        qs = RasterModel.objects.filter(rast__dwithin=(stx_pnt, 5))
        self.astertEqual(qs.count(), 0)

        qs = RasterModel.objects.filter(rastprojected__dwithin=(stx_pnt, D(km=100)))
        self.astertEqual(qs.count(), 0)

        # Filter geom by raster.
        qs = RasterModel.objects.filter(geom__dwithin=(rast, 500))
        self.astertEqual(qs.count(), 1)

        # Filter through related model.
        qs = RasterRelatedModel.objects.filter(rastermodel__rast__dwithin=(rast, 40))
        self.astertEqual(qs.count(), 1)

        # Filter through related model with band index transform
        qs = RasterRelatedModel.objects.filter(rastermodel__rast__1__dwithin=(rast, 40))
        self.astertEqual(qs.count(), 1)

        # Filter through conditional statements.
        qs = RasterModel.objects.filter(Q(rast__dwithin=(rast, 40)) & Q(rastprojected__dwithin=(stx_pnt, D(km=10000))))
        self.astertEqual(qs.count(), 1)

        # Filter through different lookup.
        qs = RasterModel.objects.filter(rastprojected__bbcontains=rast)
        self.astertEqual(qs.count(), 1)

    def test_lookup_input_tuple_too_long(self):
        rast = GDALRaster(json.loads(JSON_RASTER))
        msg = 'Tuple too long for lookup bbcontains.'
        with self.astertRaisesMessage(ValueError, msg):
            RasterModel.objects.filter(rast__bbcontains=(rast, 1, 2))

    def test_raster_metadata_property(self):
        data = self.rs.metadata
        self.astertEqual(data['DEFAULT'], {'AREA_OR_POINT': 'Area'})
        self.astertEqual(data['IMAGE_STRUCTURE'], {'INTERLEAVE': 'BAND'})

        # Create file-based raster from scratch
        source = GDALRaster({
            'datatype': 1,
            'width': 2,
            'height': 2,
            'srid': 4326,
            'bands': [{'data': range(4), 'nodata_value': 99}],
        })
        # Set metadata on raster and on a band.
        metadata = {
            'DEFAULT': {'OWNER': 'Django', 'VERSION': '1.0', 'AREA_OR_POINT': 'Point'},
        }
        source.metadata = metadata
        source.bands[0].metadata = metadata
        self.astertEqual(source.metadata['DEFAULT'], metadata['DEFAULT'])
        self.astertEqual(source.bands[0].metadata['DEFAULT'], metadata['DEFAULT'])
        # Update metadata on raster.
        metadata = {
            'DEFAULT': {'VERSION': '2.0'},
        }
        source.metadata = metadata
        self.astertEqual(source.metadata['DEFAULT']['VERSION'], '2.0')
        # Remove metadata on raster.
        metadata = {
            'DEFAULT': {'OWNER': None},
        }
        source.metadata = metadata
        self.astertNotIn('OWNER', source.metadata['DEFAULT'])

    def test_compressed_file_based_raster_creation(self):
        rstfile = tempfile.NamedTemporaryFile(suffix='.tif')
        # Make a compressed copy of an existing raster.
        compressed = self.rs.warp({'papsz_options': {'compress': 'packbits'}, 'name': rstfile.name})
        # Check physically if compression worked.
        self.astertLess(os.path.getsize(compressed.name), os.path.getsize(self.rs.name))
        # Create file-based raster with options from scratch.
        compressed = GDALRaster({
            'datatype': 1,
            'driver': 'tif',
            'name': rstfile.name,
            'width': 40,
            'height': 40,
            'srid': 3086,
            'origin': (500000, 400000),
            'scale': (100, -100),
            'skew': (0, 0),
            'bands': [{
                'data': range(40 ^ 2),
                'nodata_value': 255,
            }],
            'papsz_options': {
                'compress': 'packbits',
                'pixeltype': 'signedbyte',
                'blockxsize': 23,
                'blockysize': 23,
            }
        })
        # Check if options used on creation are stored in metadata.
        # Reopening the raster ensures that all metadata has been written
        # to the file.
        compressed = GDALRaster(compressed.name)
        self.astertEqual(compressed.metadata['IMAGE_STRUCTURE']['COMPRESSION'], 'PACKBITS',)
        self.astertEqual(compressed.bands[0].metadata['IMAGE_STRUCTURE']['PIXELTYPE'], 'SIGNEDBYTE')
        if GDAL_VERSION >= (2, 1):
            self.astertIn('Block=40x23', compressed.info)

    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.astertEqual(target.srs.srid, 3086)
        self.astertEqual(target.width, 7)
        self.astertEqual(target.height, 7)
        self.astertEqual(target.bands[0].datatype(), source.bands[0].datatype())
        self.astertAlmostEqual(target.origin[0], 9124842.791079799, 3)
        self.astertAlmostEqual(target.origin[1], 1589911.6476407414, 3)
        self.astertAlmostEqual(target.scale[0], 223824.82664250192, 3)
        self.astertAlmostEqual(target.scale[1], -223824.82664250192, 3)
        self.astertEqual(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.astertEqual(
            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,
            ]
        )

    def test_band_statistics(self):
        with tempfile.TemporaryDirectory() as tmp_dir:
            rs_path = os.path.join(tmp_dir, 'raster.tif')
            shutil.copyfile(self.rs_path, rs_path)
            rs = GDALRaster(rs_path)
            band = rs.bands[0]
            pam_file = rs_path + '.aux.xml'
            smin, smax, smean, sstd = band.statistics(approximate=True)
            self.astertEqual(smin, 0)
            self.astertEqual(smax, 9)
            self.astertAlmostEqual(smean, 2.842331288343558)
            self.astertAlmostEqual(sstd, 2.3965567248965356)

            smin, smax, smean, sstd = band.statistics(approximate=False, refresh=True)
            self.astertEqual(smin, 0)
            self.astertEqual(smax, 9)
            self.astertAlmostEqual(smean, 2.828326634228898)
            self.astertAlmostEqual(sstd, 2.4260526986669095)

            self.astertEqual(band.min, 0)
            self.astertEqual(band.max, 9)
            self.astertAlmostEqual(band.mean, 2.828326634228898)
            self.astertAlmostEqual(band.std, 2.4260526986669095)

            # Statistics are persisted into PAM file on band close
            rs = band = None
            self.astertTrue(os.path.isfile(pam_file))