From 5e13c5e385af29ae2dc570f862317ae8a038737b Mon Sep 17 00:00:00 2001 From: Jon Hill Date: Mon, 19 Nov 2018 11:03:00 +1100 Subject: [PATCH 1/3] Adding code for min/max limits. Needs testing still --- hrds/hrds.py | 41 +++++++++++++++++++++++++++++++++++------ hrds/raster.py | 15 +++++++++++++-- 2 files changed, 48 insertions(+), 8 deletions(-) diff --git a/hrds/hrds.py b/hrds/hrds.py index a444024..8ed6e66 100644 --- a/hrds/hrds.py +++ b/hrds/hrds.py @@ -38,12 +38,27 @@ 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. + 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!) + value 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 +78,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 +99,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): diff --git a/hrds/raster.py b/hrds/raster.py index 7309967..999e665 100644 --- a/hrds/raster.py +++ b/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 minmax[0] is not None: + if value < minmax[0]: + value = minmax[0] + if minmax[1] is not None: + if value > minmax[1]: + value = minmax[1] + return value @@ -196,7 +206,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? From adcb476f693eb531b81c623df0b0df38820b455b Mon Sep 17 00:00:00 2001 From: Jon Hill Date: Mon, 19 Nov 2018 11:18:52 +1100 Subject: [PATCH 2/3] Testing new min/max code. Tests passing locally --- hrds/raster.py | 17 +++++++------- tests/test_hrds.py | 44 +++++++++++++++++++++++++++++++++++++ tests/test_interpolation.py | 29 ++++++++++++++++++++++++ 3 files changed, 82 insertions(+), 8 deletions(-) diff --git a/hrds/raster.py b/hrds/raster.py index cfde93b..dc64bdb 100644 --- a/hrds/raster.py +++ b/hrds/raster.py @@ -105,12 +105,12 @@ def get_val(self, point): raise CoordinateError("Coordinate out of range", point, i, j) if self.minmax is not None: - if minmax[0] is not None: - if value < minmax[0]: - value = minmax[0] - if minmax[1] is not None: - if value > minmax[1]: - value = minmax[1] + 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 @@ -141,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:" + @@ -152,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 @@ -192,7 +193,7 @@ 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""" diff --git a/tests/test_hrds.py b/tests/test_hrds.py index ae01c3f..1600825 100644 --- a/tests/test_hrds.py +++ b/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__': diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 2acc51c..854ae57 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -46,6 +46,7 @@ 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) @@ -53,6 +54,7 @@ def test_simple_interpolation(self): 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__': From 9e0df1ae057298ebf8639d0e3cd09b744c32c187 Mon Sep 17 00:00:00 2001 From: Jon Hill Date: Mon, 19 Nov 2018 11:20:57 +1100 Subject: [PATCH 3/3] Fixing lint errors --- hrds/hrds.py | 9 +++++---- hrds/raster.py | 5 +++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/hrds/hrds.py b/hrds/hrds.py index 8ed6e66..2255fa0 100644 --- a/hrds/hrds.py +++ b/hrds/hrds.py @@ -44,10 +44,11 @@ class HRDS(object): 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!) - value when returning data. This is useful for ocean simulations - where you want a minimum depth to prevent "drying". To set this, do: + 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", diff --git a/hrds/raster.py b/hrds/raster.py index dc64bdb..0e950f3 100644 --- a/hrds/raster.py +++ b/hrds/raster.py @@ -193,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.minmax) + self.interpolator = Interpolator(origin, self.dx, self.val, + self.mask, self.minmax) def get_array(self): """return the numpy array of values""" @@ -209,7 +210,7 @@ def get_val(self, x): if (self.interpolator is None): raise RasterInterpolatorError("Should call set_band() " "before calling get_val()!") - val = self.interpolator.get_val(x) + val = self.interpolator.get_val(x) return val def point_in(self, point):