Add clip operator

README.md

@@ -86,7 +86,7 @@ Yirgacheffe is work in progress, so things planned but not supported currently:
 
 * Dynamic pixel scale adjustment - all raster layers must be provided at the same pixel scale currently  *NOW IN EXPERIMENTAL TESTING, SEE BELOW*
 * A fold operation
-* CUPY support
+* CUDA/Metal support via CUPY/MLX
 * Dispatching work across multiple CPUs *NOW IN EXPERIMENTAL TESTING, SEE BELOW*
 
 
@@ -251,6 +251,14 @@ with RasterLayer.layer_from_file('test1.tif') as layer1:
         calc.save(result)
 ```
 
+### Boolean testing
+
+Testing for equality, less than, less than or equal, greater than, and greater than or equal are supported on layers, along with logical or and logical and, as per this example, where `elevation_upper` and `elevation_lower` are scalar values:
+
+```
+filtered_elevation = (min_elevation_map <= elevation_upper) & (max_elevation_map >= elevation_lower)
+```
+
 ### Power
 
 Pixel-wise raising to a constant power:
@@ -262,6 +270,35 @@ with RasterLayer.layer_from_file('test1.tif') as layer1:
         calc.save(result)
 ```
 
+### Log, Exp, Clip, etc.
+
+The following math operators common to numpy and other libraries are currently supported:
+
+* clip
+* exp
+* exp2
+* isin
+* log
+* log2
+* log10
+* maximum
+* minimum
+* nan_to_num
+
+Typically these can be invoked either on a layer as a method:
+
+```python
+calc = layer1.log10()
+```
+
+Or via the operators module, as it's sometimes nicer to do it this way when chaining together operations in a single expression:
+
+```python
+import yirgaceffe.operators as yo
+
+calc = yo.log10(layer1 / layer2)
+```
+
 
 ### Apply
 
@@ -384,6 +421,6 @@ Because of the number of tricks that Python plays under the hood this feature ne
 
 ## Thanks
 
-Thanks to discussion and feedback from the 4C team, particularly Alison Eyres, Patrick Ferris, Amelia Holcomb, and Anil Madhavapeddy.
+Thanks to discussion and feedback from my colleagues, particularly Alison Eyres, Patrick Ferris, Amelia Holcomb, and Anil Madhavapeddy.
 
 Inspired by the work of Daniele Baisero in his AoH library.

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "yirgacheffe"
-version = "0.9.2"
+version = "0.9.3"
 description = "Abstraction of gdal datasets for doing basic math operations"
 readme = "README.md"
 authors = [{ name = "Michael Dales", email = "[email protected]" }]

tests/test_operators.py

@@ -707,7 +707,7 @@ def test_where_layers() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
-def test_isin_simple() -> None:
+def test_isin_simple_method() -> None:
     data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
     result = RasterLayer.empty_raster_layer_like(layer1)
@@ -720,6 +720,19 @@ def test_isin_simple() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
+def test_isin_simple_module() -> None:
+    data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.isin(layer1, [2, 3])
+    comp.ystep = 1
+    comp.save(result)
+
+    expected = np.isin(data1, [2, 3])
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
 def test_layer_comparison_to_value() -> None:
     data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
@@ -850,7 +863,7 @@ def test_layer_greater_than_or_equal_to_layer() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
-def test_log() -> None:
+def test_log_method() -> None:
     data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
     result = RasterLayer.empty_raster_layer_like(layer1)
@@ -862,6 +875,18 @@ def test_log() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
+def test_log_module() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.log(layer1)
+    comp.save(result)
+
+    expected = np.log(data1)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
 def test_log2() -> None:
     data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
@@ -886,7 +911,7 @@ def test_log10() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
-def test_exp() -> None:
+def test_exp_method() -> None:
     data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
     result = RasterLayer.empty_raster_layer_like(layer1)
@@ -898,6 +923,18 @@ def test_exp() -> None:
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
 
+def test_exp_module() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.exp(layer1)
+    comp.save(result)
+
+    expected = np.exp(data1)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
 def test_exp2() -> None:
     data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
@@ -937,3 +974,75 @@ def test_maximum_layers() -> None:
     expected = np.maximum(data1, data2)
     actual = result.read_array(0, 0, 4, 2)
     assert (expected == actual).all()
+
+def test_clip_both_method() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = layer1.clip(3.0, 6.0)
+    comp.save(result)
+
+    expected = data1.clip(3.0, 6.0)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
+def test_clip_both_module() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.clip(layer1, 3.0, 6.0)
+    comp.save(result)
+
+    expected = np.clip(data1, 3.0, 6.0)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
+def test_clip_upper_method() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = layer1.clip(max=6.0)
+    comp.save(result)
+
+    expected = data1.clip(max=6.0)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
+def test_clip_upper_module() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.clip(layer1, max=6.0)
+    comp.save(result)
+
+    expected = np.clip(data1, a_min=None, a_max=6.0)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
+def test_clip_lower_method() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = layer1.clip(min=3.0)
+    comp.save(result)
+
+    expected = data1.clip(min=3.0)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()
+
+def test_clip_lower_module() -> None:
+    data1 = np.array([[1.0, 2.0, 3.0, 4.0], [5.0, 2.0, 7.0, 8.0]])
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    comp = LayerOperation.clip(layer1, min=3.0)
+    comp.save(result)
+
+    expected = np.clip(data1, a_min=3.0, a_max=None)
+    actual = result.read_array(0, 0, 4, 2)
+    assert (expected == actual).all()

yirgacheffe/operators.py

@@ -75,13 +75,18 @@ def _eval(self, index, step, target_window=None):
     def nan_to_num(self, nan=0, posinf=None, neginf=None):
         return LayerOperation(
             self,
-            lambda c: np.nan_to_num(c, copy=False, nan=nan, posinf=posinf, neginf=neginf)
+            np.nan_to_num,
+            copy=False,
+            nan=nan,
+            posinf=posinf,
+            neginf=neginf,
         )
 
-    def isin(self, vals):
+    def isin(self, test_elements):
         return LayerOperation(
             self,
-            lambda c: np.isin(c, vals)
+            np.isin,
+            test_elements=test_elements,
         )
 
     def log(self):
@@ -114,6 +119,18 @@ def exp2(self):
             np.exp2
         )
 
+    def clip(self, min=None, max=None): # pylint: disable=W0622
+        # In the numpy 1 API np.clip(array) used a_max, a_min arguments and array.clip() used max and min as arguments
+        # In numpy 2 they moved so that max and min worked on both, but still support a_max, and a_min on np.clip.
+        # For now I'm only going to support the newer max/min everywhere notion, but I have to internally call
+        # a_max, a_min so that yirgacheffe can work on older numpy installs.
+        return LayerOperation(
+            self,
+            np.clip,
+            a_min=min,
+            a_max=max,
+        )
+
     def numpy_apply(self, func, other=None):
         return LayerOperation(self, func, other)
 
@@ -163,8 +180,9 @@ def minimum(a, b):
             rhs=b,
         )
 
-    def __init__(self, lhs, operator=None, rhs=None, other=None):
+    def __init__(self, lhs, operator=None, rhs=None, other=None, **kwargs):
         self.ystep = constants.YSTEP
+        self.kwargs = kwargs
 
         if lhs is None:
             raise ValueError("LHS on operation should not be none")
@@ -243,13 +261,13 @@ def _eval(self, index, step): # pylint: disable=W0221
             assert self.rhs is not None
             rhs_data = self.rhs._eval(index, step)
             other_data = self.other._eval(index, step)
-            return self.operator(lhs_data, rhs_data, other_data)
+            return self.operator(lhs_data, rhs_data, other_data, **self.kwargs)
 
         if self.rhs is not None:
             rhs_data = self.rhs._eval(index, step)
-            return self.operator(lhs_data, rhs_data)
+            return self.operator(lhs_data, rhs_data, **self.kwargs)
 
-        return self.operator(lhs_data)
+        return self.operator(lhs_data, **self.kwargs)
 
     def sum(self):
         # The result accumulator is float64, and for precision reasons
@@ -498,9 +516,9 @@ def _eval(self, index, step):
         # be nicer to promote them to arrays sooner?
         if isinstance(lhs_data, (int, float)):
             if rhs_data is None:
-                return self.operator(lhs_data)
+                return self.operator(lhs_data, **self.kwargs)
             if isinstance(rhs_data, (int, float)):
-                return self.operator(lhs_data, rhs_data)
+                return self.operator(lhs_data, rhs_data, **self.kwargs)
             else:
                 result = np.empty_like(rhs_data)
         else:
@@ -518,12 +536,21 @@ def _eval(self, index, step):
                         rhs_val = rhs_data[yoffset][xoffset]
                     except TypeError:
                         rhs_val = rhs_data
-                    result[yoffset][xoffset] = self.operator(lhs_val, rhs_val)
+                    result[yoffset][xoffset] = self.operator(lhs_val, rhs_val, **self.kwargs)
                 else:
-                    result[yoffset][xoffset] = self.operator(lhs_val)
+                    result[yoffset][xoffset] = self.operator(lhs_val, **self.kwargs)
 
         return result
 
+# We provide these module level accessors as it's often nicer to write `log(x/y)` rather than `(x/y).log()`
 where = LayerOperation.where
 minumum = LayerOperation.minimum
 maximum = LayerOperation.maximum
+clip = LayerOperation.clip
+log = LayerOperation.log
+log2 = LayerOperation.log2
+log10 = LayerOperation.log10
+exp = LayerOperation.exp
+exp2 = LayerOperation.exp2
+nan_to_num = LayerOperation.nan_to_num
+isin = LayerOperation.isin