revive ee.Image.histogramMatch

geetools/ee_image.py

@@ -1,7 +1,7 @@
 """Toolbox for the :py:class:`ee.Image` class."""
 from __future__ import annotations
 
-from typing import Any, Optional
+from typing import Any, Dict, Optional
 
 import ee
 
@@ -1246,17 +1246,34 @@ def tasseledCap(self) -> ee.Image:
     def matchHistogram(
         self,
         target: ee.Image,
-        bands: dict[str, str],
+        bands: Dict[str, str],
         geometry: ee.Geometry | None = None,
         maxBuckets: int = 256,
+        scale: int = 30,
+        maxPixels: int = 65536 * 4 - 1,
+        bestEffort: bool = True,
     ) -> ee.Image:
         """Adjust the image's histogram to match a target image.
 
+        This method performs histogram matching on individual bands independently.
+        While this adjusts the value distribution of each band to match the target,
+        it does not account for the relationships between the bands. As a result,
+        you might observe changes in the hue of the image after applying this function,
+        as the relative intensities of the different color bands are altered independently.
+
+        For histogram matching that preserves the hue of the image, consider using the
+        `matchHistogramHSV` method available in `geetools`.
+
+        From: https://medium.com/google-earth/histogram-matching-c7153c85066d
+
         Parameters:
             target: Image to match.
             bands: A dictionary of band names to match, with source bands as keys and target bands as values.
             geometry: The region to match histograms in that overlaps both images. If none is provided, the geometry of the source image will be used.
             maxBuckets: The maximum number of buckets to use when building histograms. Will be rounded to the nearest power of 2.
+            scale: the spatial resolution to compute the histograms.
+            maxPixels: The maximum number of pixels to compute the histogram.
+            bestEffort: If the polygon would contain too many pixels at the given scale, compute and use a larger scale which would allow the operation to succeed.
 
         Returns:
             The adjusted image containing the matched source bands.
@@ -1278,18 +1295,54 @@ def matchHistogram(
                 }
                 matched = source.geetools.matchHistogram(target, bands)
         """
-        raise NotImplementedError(
-            "The ee_extra package is lacking maintainer for several years, it is now incompatible with "
-            "all the latest version of Python due to use of deprecated pkg_resources. "
-            " We will try to fix this in the future, but for now please use the ee_extra package directly."
+
+        def lookup(sourceHist, targetHist):
+            """Create a lookup table to make sourceHist match targetHist."""
+            # Split the histograms by column and normalize the counts.
+            sourceValues = sourceHist.slice(1, 0, 1).project([0])
+            sourceCounts = sourceHist.slice(1, 1, 2).project([0])
+            sourceCounts = sourceCounts.divide(sourceCounts.get([-1]))  # divide each by the max
+            targetValues = targetHist.slice(1, 0, 1).project([0])
+            targetCounts = targetHist.slice(1, 1, 2).project([0])
+            targetCounts = targetCounts.divide(targetCounts.get([-1]))
+
+            # Find first position in target where targetCount >= srcCount[i], for each i.
+            lookup = sourceCounts.toList().map(
+                lambda n: targetValues.get(targetCounts.gte(n).argmax())
+            )
+            return ee.Dictionary({"x": sourceValues.toList(), "y": lookup})
+
+        bandsee = ee.Dictionary(bands)
+        geom = geometry or self._obj.geometry()
+        reducer = ee.Reducer.autoHistogram(maxBuckets=maxBuckets, cumulative=True)
+        args = dict(
+            reducer=reducer, geometry=geom, scale=scale, maxPixels=maxPixels, bestEffort=bestEffort
         )
-        # return ee_extra.Spectral.core.matchHistogram(
-        #     source=self._obj,
-        #     target=target,
-        #     bands=bands,
-        #     geometry=geometry,
-        #     maxBuckets=maxBuckets,
-        # )
+
+        # Only use pixels in target that have a value in source
+        # (inside the footprint and unmasked).
+        sourceObj = self._obj.reduceRegion(**args)
+        targetObj = target.updateMask(self._obj.mask()).reduceRegion(**args)
+
+        sourceBands = bandsee.keys()
+
+        def interpolate(sband):
+            """Interpolate bands."""
+            tband = ee.String(bandsee.get(sband))
+            lk = lookup(sourceObj.getArray(sband), targetObj.getArray(tband))
+            x = ee.List(lk.get("x"))
+            y = ee.List(lk.get("y"))
+            result = ee.Image(
+                ee.Algorithms.If(
+                    x.size().eq(1),
+                    self._obj.select([sband]).remap(x, y),
+                    self._obj.select([sband]).interpolate(x, y),
+                )
+            )
+            return result
+
+        images = sourceBands.map(interpolate)
+        return self.fromList(images)
 
     def maskClouds(
         self,

tests/conftest.py

@@ -63,6 +63,12 @@ def amazonas() -> ee.FeatureCollection:
     return colombia.filter(ee.Filter.eq("ADM1_NAME", "Amazonas"))
 
 
+@pytest.fixture
+def amazonas_centroid(amazonas) -> ee.Geometry:
+    """Return the centroid of the Amazonas state from Colombia."""
+    return ee.Feature(amazonas.first()).geometry().centroid()
+
+
 @pytest.fixture
 def s2_sr(amazonas) -> ee.ImageCollection:
     """Return a copernicus based collection.

tests/test_Image.py

@@ -285,43 +285,57 @@ def image_instance(self):
         return ee.Image(src).select(["B1", "B2", "B3"])
 
 
-class TestmatchHistogram:
+class TestMatchHistogram:
     """Test the ``histogramMatch`` method."""
 
-    @pytest.mark.xfail(reason="ee_extra package is not compatible with modern python anymore")
-    def test_histogram_match(self, image_source, image_target, vatican_buffer, num_regression):
-        bands = {"R": "R", "G": "G", "B": "B"}
+    def test_histogram_match(self, image_source, image_target, ee_image_regression):
+        """Test histogramMatch."""
+        bands = {"B8": "SR_B5", "B11": "SR_B6", "B4": "SR_B4"}
         image = image_source.geetools.matchHistogram(image_target, bands)
-        values = image.reduceRegion(ee.Reducer.mean(), vatican_buffer, 10)
-        num_regression.check(values.getInfo())
+        # for viz
+        image_target = image_target.select(["SR_B5", "SR_B6", "SR_B4"], ["N", "S", "R"])
+        image = image.select(["B8", "B11", "B4"], ["N", "S", "R"])
+        ee_image_regression.check(
+            image_target.blend(image),
+            viz_params={"bands": ["N", "S", "R"], "min": 0, "max": 0.4},
+            scale=30,
+        )
 
     @pytest.fixture
     def dates(self):
         """The dates of my imagery."""
-        return "2023-06-01", "2023-06-30"
+        return "2023-01-01", "2024-01-01"
 
     @pytest.fixture
-    def image_source(self, vatican_buffer, dates):
+    def image_source(self, amazonas_centroid, dates):
         """image from the S2 copernicus program over vatican city."""
         return (
             ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
-            .filterBounds(vatican_buffer)
+            .filter(ee.Filter.lt("CLOUD_COVERAGE_ASSESSMENT", 5))
+            .filterBounds(amazonas_centroid)
             .filterDate(*dates)
             .first()
-            .select("B4", "B3", "B2")
-            .rename("R", "G", "B")
+            .select("B8", "B11", "B4")
+            .divide(10000)
+            .toFloat()
+            .clip(amazonas_centroid.buffer(10000))
         )
 
     @pytest.fixture
-    def image_target(self, vatican_buffer, dates):
+    def image_target(self, amazonas_centroid, dates):
         """image from the L8 Landsat program over vatican city."""
         return (
             ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
-            .filterBounds(vatican_buffer)
+            .filterBounds(amazonas_centroid)
+            .filter(ee.Filter.lte("CLOUD_COVER", 5))
             .filterDate(*dates)
             .first()
-            .select("SR_B4", "SR_B3", "SR_B2")
-            .rename("R", "G", "B")
+            .select("SR_B5", "SR_B6", "SR_B4")
+            .multiply(0.0000275)
+            .add(-0.2)
+            .toFloat()
+            # .divide(65455).multiply(256).toUint8()
+            .clip(amazonas_centroid.buffer(20000).bounds())
         )