From ea6f5aa94dd8d9ff8fab9e13d675905c787d78c6 Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Fri, 27 Jan 2017 11:43:25 +0000 Subject: [PATCH] Created utility to extract pixel values from rasters by coordinates. --- docs/ref/index.rst | 1 + docs/ref/utils.rst | 35 +++++++++++++++++++++++++++ raster/utils.py | 43 +++++++++++++++++++++++++++++++++ tests/test_utils.py | 59 +++++++++++++++++++++++++++++++++++++++++++-- 4 files changed, 136 insertions(+), 2 deletions(-) create mode 100644 docs/ref/utils.rst diff --git a/docs/ref/index.rst b/docs/ref/index.rst index 8b6988e5..eda227a1 100644 --- a/docs/ref/index.rst +++ b/docs/ref/index.rst @@ -6,3 +6,4 @@ API Reference :maxdepth: 2 shortcuts + utils diff --git a/docs/ref/utils.rst b/docs/ref/utils.rst new file mode 100644 index 00000000..c967e575 --- /dev/null +++ b/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 diff --git a/raster/utils.py b/raster/utils.py index 39d18114..ac014131 100644 --- a/raster/utils.py +++ b/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] diff --git a/tests/test_utils.py b/tests/test_utils.py index 8fc98758..78da5331 100644 --- a/tests/test_utils.py +++ b/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)