Merge pull request #66 from nens/casper-reduction

Add reduction.Max

CHANGES.rst

@@ -4,6 +4,10 @@ Changelog of dask-geomodeling
 2.3.5 (unreleased)
 ------------------
 
+- Added raster.Max block to select the maximum value from overlapping sources.
+
+- Made Elemwise blocks more strict regarding source periods.
+
 - Drop python 3.5 support and move on other version requirements.
 
 - Fix deprecation warning with distutils.

dask_geomodeling/raster/__init__.py

@@ -6,3 +6,4 @@
 from .misc import *  # NOQA
 from .sources import *  # NOQA
 from .parallelize import *  # NOQA
+from .reduction import *  # NOQA

dask_geomodeling/raster/combine.py

@@ -5,17 +5,14 @@
 from datetime import timedelta as Timedelta
 import numpy as np
 
-from dask_geomodeling.utils import get_dtype_max, get_index, GeoTransform
+from dask_geomodeling.utils import filter_none, get_dtype_max, get_index
+from dask_geomodeling.utils import GeoTransform
 
 from .base import RasterBlock
 
 __all__ = ["Group"]
 
 
-def filter_none(lst):
-    return [x for x in lst if x is not None]
-
-
 class BaseCombine(RasterBlock):
     """ Base block that combines rasters into a larger one.
 

dask_geomodeling/raster/elemwise.py

@@ -55,17 +55,26 @@ def _sources(self):
         return [arg for arg in self.args if isinstance(arg, RasterBlock)]
 
     def get_sources_and_requests(self, **request):
+
+        period = self.period
+        process_kwargs = {
+            "dtype": self.dtype.name, "fillvalue": self.fillvalue,
+        }
+        if period is None:
+            return [(process_kwargs, None)]
+
+        # limit request to self.period so that resulting data is aligned
         start = request.get("start", None)
         stop = request.get("stop", None)
-
-        if start is not None and stop is not None:
-            # limit request to self.period so that resulting data is aligned
-            period = self.period
-            if period is not None:
+        if start is not None:
+            if stop is not None:
                 request["start"] = max(start, period[0])
                 request["stop"] = min(stop, period[1])
+            else:  # stop is None but start isn't
+                request["start"] = min(max(start, period[0]), period[1])
+        else:  # start and stop are both None
+            request["start"] = period[1]
 
-        process_kwargs = {"dtype": self.dtype.name, "fillvalue": self.fillvalue}
         sources_and_requests = [(source, request) for source in self.args]
 
         return [(process_kwargs, None)] + sources_and_requests
@@ -241,6 +250,9 @@ def wrap_math_process_func(func):
 
     @wraps(func)  # propagates the name and docstring
     def math_process_func(process_kwargs, *args):
+        if not args:
+            return None
+
         compute_args = []  # the args for the math operation
         # perform the nodata masking manually, as numpy maskedarrays are slow
         nodata_mask = None

dask_geomodeling/raster/reduction.py

@@ -2,10 +2,14 @@
 Module containing reduction raster blocks.
 """
 import numpy as np
-from dask_geomodeling.utils import get_index, parse_percentile_statistic
-
+from dask_geomodeling.utils import filter_none, get_index
+from dask_geomodeling.utils import parse_percentile_statistic
+from .base import RasterBlock
+from .elemwise import BaseElementwise
 from functools import partial
 
+__all__ = ["Max"]
+
 STATISTICS = {
     "first": None,
     "last": None,
@@ -115,3 +119,80 @@ def reduce_rasters(stack, statistic, no_data_value=None, dtype=None):
         out[not_all_nan] = func(stack_array[:, not_all_nan], axis=0)
 
     return {"values": out, "no_data_value": no_data_value}
+
+
+class BaseReduction(BaseElementwise):
+    def __init__(self, *args):
+        for arg in args:
+            if not isinstance(arg, RasterBlock):
+                raise TypeError("'{}' object is not allowed".format(type(arg)))
+        super().__init__(*args)
+
+    @property
+    def extent(self):
+        """ Boundingbox of combined contents in WGS84 projection. """
+        extents = filter_none([x.extent for x in self.args])
+        if len(extents) == 0:
+            return None
+        elif len(extents) == 1:
+            return extents[0]
+
+        # multiple extents: return the joined box
+        x1 = min([e[0] for e in extents])
+        y1 = min([e[1] for e in extents])
+        x2 = max([e[2] for e in extents])
+        y2 = max([e[3] for e in extents])
+        return x1, y1, x2, y2
+
+    @property
+    def geometry(self):
+        """Combined geometry in the projection of the first store geometry. """
+        geometries = filter_none([x.geometry for x in self.args])
+        if len(geometries) == 0:
+            return
+        elif len(geometries) == 1:
+            return geometries[0]
+        result = geometries[0]
+        sr = result.GetSpatialReference()
+        for geometry in geometries[1:]:
+            if not geometry.GetSpatialReference().IsSame(sr):
+                geometry = geometry.Clone()
+                geometry.TransformTo(sr)
+            result = result.Union(geometry)
+        return result
+
+
+def wrap_reduction_function(statistic):
+    def reduction_function(process_kwargs, *args):
+        # remove None values
+        stack = [x for x in args if x is not None]
+        # return None if all source data is None
+        if len(stack) == 0:
+            return
+
+        # see BaseElementWise.get_sources_and_requests
+        return reduce_rasters(
+            stack,
+            statistic,
+            process_kwargs["fillvalue"],
+            process_kwargs["dtype"],
+        )
+    return reduction_function
+
+
+class Max(BaseReduction):
+    """
+    Take the maximum value of two or more rasters, ignoring no data.
+
+    Args:
+      *args (list of RasterBlocks): list of rasters to be combined.
+
+    Returns:
+      RasterBlock with the maximum values
+    """
+    process = staticmethod(wrap_reduction_function("max"))
+
+    @property
+    def dtype(self):
+        # skip the default behaviour where we use at least int32 / float32
+        return np.result_type(*self.args)

dask_geomodeling/tests/test_raster_reduction.py

@@ -3,6 +3,7 @@
 import pytest
 
 from dask_geomodeling.raster.reduction import reduce_rasters
+from dask_geomodeling import raster
 
 
 @pytest.fixture
@@ -96,3 +97,21 @@ def test_reduce_defaults(statistic, stack):
 def test_reduce_raises_zero_length(stack):
     with pytest.raises(ValueError):
         reduce_rasters([], "first")
+
+
+def test_max(source, vals_request):
+    block = raster.Max(source, source - 2)
+    data = block.get_data(**vals_request)
+    assert data["values"][:, 0, 0].tolist() == [1, 7, data["no_data_value"]]
+
+
+def test_max_with_nodata(source, nodata_source, vals_request):
+    block = raster.Max(source, nodata_source)
+    data = block.get_data(**vals_request)
+    assert data["values"][:, 0, 0].tolist() == [1, 7, data["no_data_value"]]
+
+
+def test_max_with_empty(source, empty_source, vals_request):
+    block = raster.Max(source, empty_source)
+    data = block.get_data(**vals_request)
+    assert data is None

dask_geomodeling/utils.py

@@ -941,3 +941,7 @@ def dt_to_ms(dt):
     if dt.tzinfo is None:
         dt = dt.replace(tzinfo=pytz.UTC)
     return int(dt.timestamp() * 1000)
+
+
+def filter_none(lst):
+    return [x for x in lst if x is not None]