Created utility to extract pixel values from rasters by coordinates.

docs/ref/index.rst

@@ -6,3 +6,4 @@ API Reference
    :maxdepth: 2
 
    shortcuts
+   utils

docs/ref/utils.rst

@@ -0,0 +1,35 @@
+================
+Raster Utilities
+================
+Django-raster hosts some utilities that ease the interaction with raster data.
+The functions are located in ``raster.utils`` and ``raster.tiles.utils``.
+
+
+.. function:: pixel_value_from_point(raster, point, band=0)
+
+    Return the pixel value for the coordinate of the input point from selected
+    band. 
+
+    The input can be a point or tuple, if its a tuple it is assumed to be a
+    pair of coordinates in the reference system of the raster. The band index
+    to be used for extraction can be specified with the ``band`` keyword.
+
+    Example::
+
+        # Create a raster.
+        >>> raster = GDALRaster({
+            'width': 5,
+            'height': 5,
+            'srid': 4326,
+            'bands': [{'data': range(25)}],
+            'origin': (2, 2),
+            'scale': (1, 1)
+        })
+        # Create a point at origin
+        >>> point = OGRGeometry('SRID=4326;POINT(2 2)')
+        # Get pixel value at origin.
+        >>> pixel_value_from_point(raster, point)
+        ... 0
+        # Get pixel value from within the raster, using coordinate tuple input.
+        >>> pixel_value_from_point(raster, (2, 3.5))
+        ... 5

raster/utils.py

@@ -3,6 +3,7 @@
 import numpy
 from PIL import Image
 
+from django.contrib.gis.gdal import OGRGeometry
 from django.utils import six
 from raster.algebra.parser import FormulaParser
 from raster.exceptions import RasterException
@@ -129,3 +130,45 @@ def colormap_to_rgba(colormap):
         return {k: hex_to_rgba(v) if isinstance(v, (six.string_types, int)) and k in ['from', 'to', 'over'] else v for k, v in colormap.items()}
     else:
         return {k: hex_to_rgba(v) if isinstance(v, (six.string_types, int)) else v for k, v in colormap.items()}
+
+
+def pixel_value_from_point(raster, point, band=0):
+    """
+    Returns the pixel value for the coordinate of the input point from selected
+    band.
+
+    The input can be a point or tuple, if its a tuple it is assumed to be
+    coordinates in the reference system of the raster.
+    """
+    if isinstance(point, (tuple, list)):
+        point = OGRGeometry('POINT({0} {1})'.format(*point))
+        point.srid = raster.srid
+    elif not point.srs or not raster.srs:
+        raise ValueError('Both the point and the raster are required to have a reference system specified.')
+    elif point.srs != raster.srs:
+        # Ensure the projection of the point is the same as of the raster.
+        point.transform(raster.srid)
+
+    # Return if point and raster do not touch.
+    bbox = OGRGeometry.from_bbox(raster.extent)
+    bbox.srs = raster.srs
+
+    if not point.intersects(bbox):
+        return
+
+    # Compute position of point relative to raster origin.
+    offset = (abs(raster.origin.x - point.coords[0]), abs(raster.origin.y - point.coords[1]))
+
+    # Compute pixel index value based on offset.
+    offset_index = [int(offset[0] / abs(raster.scale.x)), int(offset[1] / abs(raster.scale.y))]
+
+    # If the point is exactly on the boundary, the offset_index is rounded to
+    # a pixel index over the edge of the pixel. The index needs to be reduced
+    # by one pixel for those cases.
+    if offset_index[0] == raster.width:
+        offset_index[0] -= 1
+
+    if offset_index[1] == raster.height:
+        offset_index[1] -= 1
+
+    return raster.bands[band].data(offset=offset_index, size=(1, 1))[0, 0]

tests/test_utils.py

@@ -2,11 +2,11 @@
 
 import numpy
 
-from django.contrib.gis.gdal import OGRGeometry
+from django.contrib.gis.gdal import GDALRaster, OGRGeometry
 from django.test import TestCase
 from raster.exceptions import RasterException
 from raster.tiles.utils import tile_bounds, tile_index_range
-from raster.utils import colormap_to_rgba, hex_to_rgba, rescale_to_channel_range
+from raster.utils import colormap_to_rgba, hex_to_rgba, pixel_value_from_point, rescale_to_channel_range
 
 
 class TestUtils(TestCase):
@@ -88,3 +88,58 @@ def test_channel_rescale(self):
             rescale_to_channel_range(data, 60, 50, 40),
             [60, 40, 50]
         )
+
+    def test_get_pixel_value(self):
+        raster = GDALRaster({'width': 5, 'height': 5, 'srid': 4326, 'bands': [{'data': range(25)}], 'origin': (2, 2), 'scale': (1, 1)})
+
+        # Pixel value at origin.
+        point = OGRGeometry('SRID=4326;POINT(2 2)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 0)
+
+        # Coords as tuple.
+        result = pixel_value_from_point(raster, (2, 2))
+        self.assertEqual(result, 0)
+
+        # Point in different projection.
+        point = OGRGeometry('SRID=3857;POINT(222638.9815865472 222684.20850554455)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 0)
+
+        # Pixel value outside of raster.
+        point = OGRGeometry('SRID=4326;POINT(-2 2)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, None)
+
+        point = OGRGeometry('SRID=4326;POINT(8 8)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, None)
+
+        # Pixel values within the raster.
+        point = OGRGeometry('SRID=4326;POINT(3.5 2)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 1)
+
+        point = OGRGeometry('SRID=4326;POINT(2 3.5)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 5)
+
+        point = OGRGeometry('SRID=4326;POINT(6.999 6.9999)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 24)
+
+        # Pixel value at "outer" edge of raster.
+        point = OGRGeometry('SRID=4326;POINT(7 7)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 24)
+
+        # Point without srs specified.
+        point = OGRGeometry('POINT(2 2)')
+        with self.assertRaises(ValueError):
+            pixel_value_from_point(raster, point)
+
+        # Raster with negative scale on y axis.
+        raster = GDALRaster({'width': 5, 'height': 5, 'srid': 4326, 'bands': [{'data': range(25)}], 'origin': (2, 2), 'scale': (1, -1)})
+        point = OGRGeometry('SRID=4326;POINT(3 1)')
+        result = pixel_value_from_point(raster, point)
+        self.assertEqual(result, 6)