Merge pull request #6 from EnvModellingGroup/nan_distance

Add ability to create distance based on missing data (nodata) in the rasters, as well as the edge. This allows for partial spatial data sets, e.g. processed bathymetry surveys, rather than just DEM "products", such as GEBCO. Fixes #4

hrds/raster.py

@@ -141,6 +141,7 @@ def __init__(self, filename):
         self.interpolator = None
         self.extent = None
         self.dx = 0.0
+        self.nodata = None
 
     def get_extent(self):
         """Return list of corner coordinates from a geotransform
@@ -175,6 +176,7 @@ def set_band(self, band_no=1):
            Usually 1, which is default"""
         self.band = band_no
         raster = self.ds.GetRasterBand(self.band)
+        self.nodata = raster.GetNoDataValue()
         self.val = np.flipud(np.array(raster.ReadAsArray()))
         self.extent = self.get_extent()
         origin = np.amin(self.extent, axis=0)

hrds/raster_buffer.py

@@ -105,8 +105,15 @@ def make_buffer(self, output_file):
 
         # fill with edge value
         dist = np.full((ncols, nrows), 0.0)
-        # then fill in the middle
+        # then fill in the middle,
+        # except where there is no data
         dist[1:-1, 1:-1] = 1
+        if self.over is None:
+            orig_raster = self.raster.get_array()
+            nodata = self.raster.nodata
+            # TODO this will only work if dist and orig_raster
+            # are the same size
+            dist[orig_raster == nodata] = 0
         # calc euclidian distance and convert to 0 -> 1 scale
         dist = distance_transform_edt(dist, sampling=[dx[0], dx[1]])
         dist = dist / self.distance

tests/test_buffer.py

@@ -19,6 +19,7 @@
 # Copyright Jon Hill, University of York, [email protected]
 
 test_file_name1 = "tests/test_raster_large.tif"
+test_file_name2 = "tests/test_raster_nodata_centre.tif"
 temp_file = "temp.tif"
 
 class TestBufferCreator(unittest.TestCase):
@@ -87,5 +88,32 @@ def test_simple_distance_setres(self):
         self.assertEqual(rci.get_val(point3),1.0)
         self.assertAlmostEqual(rci.get_val(point4),0.5,delta=0.01)
 
+    def test_simple_distance_nan(self):
+        """ Very simple test with a raster object like thus:
+             1  2  3  4
+             5  6  7  8
+             9  NAN NAN 12
+             13 14 15 16
+            (on a 20x20 grid, not 4x4 as above...only so much space 
+            in comments!). The slightly-left-of-centre has been 
+            replaced by NaN
+            
+            LLC is 0,0 and upper right is 4,4. (hence dx is 0.2)
+            The data are stored in cell centres and we ask for a few coords.
+            We check that the value is 1 at the distance, 0 at the boundary
+            and, for this test, the value is 0 in the centre.
+            """
+        rbuff = CreateBuffer(test_file_name2, 1.0)
+        point1 = [0.2, 3.85] # should return ~0
+        point2 = [0.7, 2] # should be 0.4
+        point3 = [1.7, 2] # should be 0 (in the centre)
+        rbuff.make_buffer(temp_file)
+        # we now read in the buffer using the rasterinterpolator class
+        rci = RasterInterpolator(temp_file)
+        rci.set_band()
+        self.assertAlmostEqual(rci.get_val(point1),0.0, delta=0.03)
+        self.assertAlmostEqual(rci.get_val(point2),0.4, delta=0.03)
+        self.assertEqual(rci.get_val(point3),0.0)
+
 if __name__ == '__main__':
     unittest.main()

tests/test_hrds.py

@@ -54,7 +54,8 @@ def test_simple_rasters(self):
         for p,e in zip(points, expected):
             self.assertAlmostEqual(bathy.get_val(p),e,delta=0.1)
 
-@unittest.skip("Skipping this by default. Uses proprietary data.")
+@unittest.skipUnless(os.path.isfile("tests/real_data/gebco_uk.tif"),
+                 "Skipping as proprietary data missing.")
 class RealDataTest(unittest.TestCase):
 
     def test_real_world(self):

tests/test_interpolation.py

@@ -78,7 +78,8 @@ def test_point_in(self):
         self.assertTrue(rci.point_in(point4))
         self.assertFalse(rci.point_in(point5))
 
-    @unittest.skip("Skipping this by default. Uses proprietary data.")
+    @unittest.skipUnless(os.path.isfile("tests/real_data/gebco_uk.tif"),
+                         "Skipping as proprietary data missing.")
     def test_real_data(self):
         """Using gebco cut out to test interpolation
         """