Merge pull request #71 from nens/verkerk-logexp

Add Exp and Log RasterBlocks

CHANGES.rst

@@ -4,7 +4,7 @@ Changelog of dask-geomodeling
 2.2.13 (unreleased)
 -------------------
 
-- Nothing changed yet.
+- Added Exp, Log and Log10 RasterBlocks.
 
 
 2.2.12 (2020-09-29)

dask_geomodeling/raster/elemwise.py

@@ -28,6 +28,9 @@
     "Xor",
     "IsData",
     "IsNoData",
+    "Exp",
+    "Log",
+    "Log10",
 ]
 
 
@@ -283,6 +286,9 @@ def math_process_func(process_kwargs, *args):
         with np.errstate(all="ignore"):  # suppresses warnings
             result_values = func(*compute_args, **func_kwargs)
 
+        # mask nan & inf values
+        result_values[~np.isfinite(result_values)] = fillvalue
+
         if nodata_mask is not None:
             result_values[nodata_mask] = fillvalue
         return {"no_data_value": no_data_value, "values": result_values}
@@ -304,7 +310,7 @@ class Add(BaseMath):
 
     Returns:
       RasterBlock containing the result of the addition.
-	"""
+    """
 
     process = staticmethod(wrap_math_process_func(np.add))
 
@@ -339,7 +345,7 @@ class Multiply(BaseMath):
     Args:
       a (RasterBlock, number): Multiplication factor a
       b (RasterBlock, number): Multiplication factor b
-     
+
     Returns:
       RasterBlock containing the result of the multiplication.
     """
@@ -358,7 +364,7 @@ class Divide(BaseMath):
     Args:
       a (RasterBlock, number): Numerator
       b (RasterBlock, number): Denominator
-     
+
     Returns:
       RasterBlock containing the result of the division.
     """
@@ -382,9 +388,9 @@ class Power(BaseMath):
     Args:
       a (RasterBlock, number): Base
       b (RasterBlock, number): Exponent
-      
+
     Returns:
-      RasterBlock containing the result of the exponential function. 
+      RasterBlock containing the result of the exponential function.
     """
 
     process = staticmethod(wrap_math_process_func(np.power))
@@ -412,7 +418,7 @@ class Equal(BaseComparison):
     Args:
       a (RasterBlock, number): Comparison term a
       b (RasterBlock, number): Comparison term b
-      
+
     Returns:
       RasterBlock containing boolean values
     """
@@ -479,7 +485,7 @@ class GreaterEqual(BaseComparison):
     Either one or both of the inputs should be a RasterBlock. In case of
     two raster inputs the temporal properties of the rasters should be equal,
     however spatial properties can be different.
-    
+
     Args:
       a (RasterBlock, number): Comparison term a
       b (RasterBlock, number): Comparison term b
@@ -503,7 +509,7 @@ class Less(BaseComparison):
     Either one or both of the inputs should be a RasterBlock. In case of
     two raster inputs the temporal properties of the rasters should be equal,
     however spatial properties can be different.
-    
+
     Args:
       a (RasterBlock, number): Comparison term a
       b (RasterBlock, number): Comparison term b
@@ -527,7 +533,7 @@ class LessEqual(BaseComparison):
     Either one or both of the inputs should be a RasterBlock. In case of
     two raster inputs the temporal properties of the rasters should be equal,
     however spatial properties can be different.
-    
+
     Args:
       a (RasterBlock, number): Comparison term a
       b (RasterBlock, number): Comparison term b
@@ -542,15 +548,15 @@ class LessEqual(BaseComparison):
 class Invert(BaseSingle):
     """
     Logically invert a raster (swap True and False).
-    
+
     Takes a single input raster containing boolean values and outputs a boolean
     raster with the same spatial and temportal properties.
 
     Args:
       x (RasterBlock): Boolean raster with values to invert
 
     Returns:
-      RasterBlock with boolean values opposite to the input raster. 
+      RasterBlock with boolean values opposite to the input raster.
     """
 
     def __init__(self, x):
@@ -573,15 +579,15 @@ def dtype(self):
 class IsData(BaseSingle):
     """
     Returns True where raster has data.
-    
+
     Takes a single input raster and outputs a boolean raster with the same
     spatial and temporal properties.
 
     Args:
       store (RasterBlock): Input raster
 
     Returns:
-      RasterBlock with boolean values. 
+      RasterBlock with boolean values.
     """
 
     def __init__(self, store):
@@ -617,7 +623,7 @@ class IsNoData(IsData):
       store (RasterBlock): Input raster
 
     Returns:
-      RasterBlock with boolean values. 
+      RasterBlock with boolean values.
     """
 
     @staticmethod
@@ -636,7 +642,7 @@ class And(BaseLogic):
     Either one or both of the inputs should be a boolean RasterBlock. In case
     of two raster inputs the temporal properties of the rasters should be
     equal, however spatial properties can be different.
-    
+
     Args:
       a (RasterBlock, boolean): Logical term a
       b (RasterBlock, boolean): Logical term b
@@ -702,7 +708,7 @@ class FillNoData(BaseElementwise):
 
     Args:
       *args (list of RasterBlocks): Rasters to be combined.
-      
+
     Returns:
       RasterBlock that combines values from the inputs.
     """
@@ -746,3 +752,74 @@ def process(process_kwargs, *args):
             values[index] = data[index]
 
         return {"values": values, "no_data_value": fillvalue}
+
+
+class BaseLogExp(BaseSingle):
+    """
+    Base Block for Exp, Log and Log10 Blocks.
+    """
+    def __init__(self, x):
+        if x.dtype == np.dtype("bool"):
+            raise TypeError("input block must not have boolean dtype")
+        super(BaseLogExp, self).__init__(x)
+
+    def get_sources_and_requests(self, **request):
+        process_kwargs = {"dtype": self.dtype.name, "fillvalue": self.fillvalue}
+        return [(process_kwargs, None), (self.args[0], request)]
+
+    @property
+    def dtype(self):
+        # use at least float32
+        return np.result_type(np.float32, *self.args)
+
+    @property
+    def fillvalue(self):
+        return get_dtype_max(self.dtype)
+
+
+class Exp(BaseLogExp):
+    """
+    Return e raised to the power of the raster values.
+
+    Out-of-range results (not representable by the resulting datatype) are set
+    to `no data`.
+
+    Args:
+      x (RasterBlock): Raster
+
+    Returns:
+      RasterBlock.
+    """
+    process = staticmethod(wrap_math_process_func(np.exp))
+
+
+class Log(BaseLogExp):
+    """
+    Return natural logarithm of the raster values.
+
+    Out-of-range results (not representable by the resulting datatype) are set
+    to `no data` as well as the result of input values < 0.
+
+    Args:
+      x (RasterBlock): Raster
+
+    Returns:
+      RasterBlock.
+    """
+    process = staticmethod(wrap_math_process_func(np.log))
+
+
+class Log10(BaseLogExp):
+    """
+    Return the base 10 logarithm of the raster values.
+
+    Out-of-range results (not representable by the resulting datatype) are set
+    to `no data` as well as the result of input values < 0.
+
+    Args:
+      x (RasterBlock): Raster
+
+    Returns:
+      RasterBlock.
+    """
+    process = staticmethod(wrap_math_process_func(np.log10))

dask_geomodeling/tests/factories.py

@@ -49,7 +49,7 @@ def dtype(self):
         try:
             return self.value.dtype
         except AttributeError:
-            return np.uint8
+            return np.dtype(np.uint8)
 
     @property
     def fillvalue(self):

dask_geomodeling/tests/test_raster.py

@@ -7,7 +7,8 @@
 from scipy import ndimage
 
 from dask_geomodeling import raster
-from dask_geomodeling.utils import EPSG4326, EPSG3857, Extent, get_epsg_or_wkt
+from dask_geomodeling.utils import EPSG4326, EPSG3857
+from dask_geomodeling.utils import Extent, get_dtype_max, get_epsg_or_wkt
 from dask_geomodeling.raster import RasterBlock
 from dask_geomodeling.tests.factories import MockRaster, MockGeometry
 
@@ -292,6 +293,15 @@ def setUp(self):
             bands=1,
             value=[[1, 1], [7, 7], [255, 255]],  # 2x3 shape. 255 is nodata
         )
+        self.logexp_storage = MockRaster(
+            origin=Datetime(2000, 1, 1),
+            timedelta=Timedelta(hours=1),
+            bands=1,
+            value=np.array([
+                [-1, 0],
+                [np.e, 10],
+                [999, get_dtype_max('f8')],
+            ], dtype='f8'))
         self.vals_request = dict(
             mode="vals",
             start=Datetime(2000, 1, 1),
@@ -538,6 +548,35 @@ def test_math_nodata(self):
             result = view.get_data(**self.vals_request)
             assert_equal(result["values"], result["no_data_value"])
 
+    def test_base_log_exp_init(self):
+        # cannot take Exp from boolean storage
+        arg = self.logexp_storage == 7
+        self.assertRaises(TypeError, raster.elemwise.BaseLogExp, arg)
+
+    def test_exp(self):
+        view = raster.Exp(self.logexp_storage)
+        n = view.fillvalue
+        expected = [[1 / np.e, 1], [np.exp(np.e), np.exp(10)], [n, n]]
+        assert_allclose(
+            view.get_data(**self.vals_request)["values"][0], expected,
+        )
+
+    def test_log_e(self):
+        view = raster.Log(self.logexp_storage)
+        n = view.fillvalue
+        expected = [[n, n], [1, np.log(10)], [np.log(999), n]]
+        assert_allclose(
+            view.get_data(**self.vals_request)["values"][0], expected,
+        )
+
+    def test_log_10(self):
+        view = raster.Log10(self.logexp_storage)
+        n = view.fillvalue
+        expected = [[n, n], [np.log10(np.e), 1], [np.log10(999), n]]
+        assert_allclose(
+            view.get_data(**self.vals_request)["values"][0], expected,
+        )
+
 
 class TestFillNoData(unittest.TestCase):
     klass = raster.FillNoData