Merge pull request #7 from EnvModellingGroup/min-max_depth

Adding capability of specifying a min/max depth per raster. Fixes #5

hrds/hrds.py

@@ -38,12 +38,28 @@ class HRDS(object):
     bathy = HRDS("gebco_uk.tif",
              rasters=("emod_utm.tif",
                       "marine_digimap.tif"),
-             distances=(10000, 5000))
+             distances=(10000, 5000),
+             minmax=None)
     bathy.set_bands()
 
     The first argument is the base raster filename. `rasters` is a list
     of raster filenames, with corresponding `distances` over which to
     create the buffer. Distances are in the same units as the rasters.
+    The min/max argument allows you to specify a minimum or maximum
+    (or both!) values when returning data. This is useful for ocean
+    simulations where you want a minimum depth to prevent "drying".
+    To set this, do:
+
+    bathy = HRDS("gebco_uk.tif",
+             rasters=("emod_utm.tif",
+                      "marine_digimap.tif"),
+             distances=(10000, 5000),
+             minmax=[[None,-5],[None,-3],[None,None]])
+
+    which would set a maximum depth of -5m on the gebco data (-ve = below
+    sea level, +ve above), maximum of -3m on the emod data and no limits
+    on the marine_digimap data. You must supply the same number of min-max
+    pairs are there are total rasters.
 
     If is possible to supply the buffer rasters directly (e.g. if you want
     to use different distances on each edge of your raster, or some other
@@ -63,7 +79,8 @@ class HRDS(object):
     bathy.get_val(100,100)
 
         """
-    def __init__(self, baseRaster, rasters=None, distances=None, buffers=None):
+    def __init__(self, baseRaster, rasters=None, distances=None,
+                 buffers=None, minmax=None):
         """ baseRaster is the low res raster filename across whole domain.
         rasters is a list of filenames of the other rasters in priority order.
         distances is the distance to create a buffer (in same units as
@@ -83,10 +100,23 @@ def __init__(self, baseRaster, rasters=None, distances=None, buffers=None):
                                 "rasters and "+str(len(distances)) +
                                 "distances. They should match")
 
-        self.baseRaster = RasterInterpolator(baseRaster)
+        if minmax is not None:
+            if len(rasters)+1 != len(minmax):
+                raise HRDSError("Please supply same number of minmax values" +
+                                "as the total number of rasters, inc. base." +
+                                "You gave me: "+str(len(minmax))+" min/max" +
+                                "and I expected: "+str(len(rasters)+1))
+
+        if minmax is None:
+            self.baseRaster = RasterInterpolator(baseRaster)
+        else:
+            self.baseRaster = RasterInterpolator(baseRaster, minmax[0])
         self.raster_stack = []
-        for r in rasters:
-            self.raster_stack.append(RasterInterpolator(r))
+        for i, r in enumerate(rasters):
+            if minmax is None:
+                self.raster_stack.append(RasterInterpolator(r))
+            else:
+                self.raster_stack.append(RasterInterpolator(r, minmax[i+1]))
         self.buffer_stack = []
         if buffers is None:
             for r, d in zip(rasters, distances):

hrds/raster.py

@@ -45,11 +45,12 @@ class Interpolator(object):
     may switch bands and hence have to reload the val data.
     """
 
-    def __init__(self, origin, delta, val, mask=None):
+    def __init__(self, origin, delta, val, mask=None, minmax=None):
         self.origin = origin
         self.delta = delta
         self.val = val
         self.mask = mask
+        self.minmax = minmax
 
     def set_mask(self, mask):
         self.mask = mask
@@ -102,6 +103,15 @@ def get_val(self, point):
                                                   "should have 2 dimensions")
         except IndexError:
             raise CoordinateError("Coordinate out of range", point, i, j)
+
+        if self.minmax is not None:
+            if self.minmax[0] is not None:
+                if value < self.minmax[0]:
+                    value = self.minmax[0]
+            if self.minmax[1] is not None:
+                if value > self.minmax[1]:
+                    value = self.minmax[1]
+
         return value
 
 
@@ -131,7 +141,7 @@ class RasterInterpolator(object):
     calls of set_band().
 
     """
-    def __init__(self, filename):
+    def __init__(self, filename, minmax=None):
         self.ds = gdal.Open(filename)
         if (self.ds is None):
             raise RasterInterpolatorError("Couldn't find your raster file:" +
@@ -142,6 +152,7 @@ def __init__(self, filename):
         self.extent = None
         self.dx = 0.0
         self.nodata = None
+        self.minmax = minmax
 
     def get_extent(self):
         """Return list of corner coordinates from a geotransform
@@ -182,7 +193,8 @@ def set_band(self, band_no=1):
         origin = np.amin(self.extent, axis=0)
         transform = self.ds.GetGeoTransform()
         self.dx = [transform[1], -transform[5]]
-        self.interpolator = Interpolator(origin, self.dx, self.val, self.mask)
+        self.interpolator = Interpolator(origin, self.dx, self.val,
+                                         self.mask, self.minmax)
 
     def get_array(self):
         """return the numpy array of values"""
@@ -198,7 +210,8 @@ def get_val(self, x):
         if (self.interpolator is None):
             raise RasterInterpolatorError("Should call set_band() "
                                           "before calling get_val()!")
-        return self.interpolator.get_val(x)
+        val = self.interpolator.get_val(x)
+        return val
 
     def point_in(self, point):
         # does this point occur in the raster?

tests/test_hrds.py

@@ -97,6 +97,50 @@ def test_real_world(self):
         for p,e in zip(points, expected):
             self.assertAlmostEqual(bathy.get_val(p),e,delta=0.75)
 
+
+@unittest.skipUnless(os.path.isfile("tests/real_data/gebco_uk.tif"),
+                 "Skipping as proprietary data missing.")
+class RealDataTest(unittest.TestCase):
+
+    def test_real_world_limit(self):
+        """ Mix of GEBCO, EMODnet and Marine Digimap just off the
+        Norfolk/Suffolk coast.
+        Projected to UTM 30, so everything in m.
+        """
+        bathy = HRDS("tests/real_data/gebco_uk.tif", 
+                     rasters=("tests/real_data/emod_utm.tif", 
+                              "tests/real_data/marine_digimap.tif"), 
+                     distances=(10000, 5000),
+                     minmax=[[None,-25],[None,-30],[None,None]])
+        bathy.set_bands()
+        """
+        Our data:
+        X          Y      emod_utm  gebco_uk  marine_dig   ID
+        823862.	5782011.           -21.21704                1
+        839323.	5782408.  -25.4705 -24.032                  2
+        853000.	5782804.  -43.1108 -38.03058                3
+        858947.	5782606.  -50.5894 -46.71868   -52.03551    4
+        866083.	5783201.  -43.4241 -40.0147    -48.12579    5
+        889870.	5784787.  -41.1196 -32.12536                6
+        949138.	5782408.           -22.1890                 7
+"""
+        points = ([823862., 5782011.],
+                  [839323., 5782408.],
+                  [853000., 5782804.],                  
+                  [858947., 5782606.],
+                  [866083., 5783201.],
+                  [889870., 5784787.],
+                  [949138., 5782408.],
+                  )
+        expected = [-25.0,  # -21.21704 - limited!
+                    -25.5, # estimate and partially limited
+                    -43.1108,
+                    -50.0, # estimate
+                    -48.12579,
+                    -41.1196,
+                    -25.0] #-22.189 - limited
+        for p,e in zip(points, expected):
+            self.assertAlmostEqual(bathy.get_val(p),e,delta=0.75)
 
 
 if __name__ == '__main__':

tests/test_interpolation.py

@@ -46,13 +46,15 @@ def test_simple_interpolation(self):
         point4 = [3,1] # should return 13.5
         point5 = [1.999999, 2.999999] # should return nearly 4.5
         point6 = [0.5, 0.5] # should return 13 (llc)
+        point7 = [1.0, 3.49] # should be 1.5
         self.assertRaises(RasterInterpolatorError, rci.get_val, point2)
         rci.set_band()        
         self.assertEqual(rci.get_val(point2),8.0)
         self.assertAlmostEqual(rci.get_val(point3),4.5)
         self.assertAlmostEqual(rci.get_val(point5),4.5,5)
         self.assertEqual(rci.get_val(point4),13.5)
         self.assertEqual(rci.get_val(point6),13)
+        self.assertAlmostEqual(rci.get_val(point7),1.5, delta=0.1)
         self.assertRaises(CoordinateError,rci.get_val, point1)
 
     def test_point_in(self):
@@ -104,6 +106,33 @@ def test_real_data(self):
         for p,e in zip(points, expected):
             self.assertAlmostEqual(rci.get_val(p),e,delta=0.75)
 
+    def test_simple_interpolation_limits(self):
+        """ Very simple test with a raster object like thus:
+             1  2  3  4
+             5  6  7  8
+             9  10 11 12
+             13 14 15 16
+            
+            LLC is 0,0 and upper right is 4,4. 
+            But here we set upper and lower limits
+            """
+        rci = RasterInterpolator(test_file_name1,minmax=[2,10])
+        point1 = [0.0, 0.0] # should error
+        point2 = [1.5, 2.0] # should return 8
+        point3 = [2.0, 3.0] # should return 4.5
+        point4 = [3,1] # should return 13.5, limited to 10
+        point5 = [1.999999, 2.999999] # should return nearly 4.5
+        point6 = [0.5, 0.5] # should return 13 (llc), limited to 10
+        point7 = [1.0, 3.49] # should be 1.5, but limited to 2
+        self.assertRaises(RasterInterpolatorError, rci.get_val, point2)
+        rci.set_band()        
+        self.assertEqual(rci.get_val(point2),8.0)
+        self.assertAlmostEqual(rci.get_val(point3),4.5)
+        self.assertAlmostEqual(rci.get_val(point5),4.5,5)
+        self.assertEqual(rci.get_val(point4),10)
+        self.assertEqual(rci.get_val(point6),10)
+        self.assertEqual(rci.get_val(point7),2)
+        self.assertRaises(CoordinateError,rci.get_val, point1)
 
 
 if __name__ == '__main__':