Add a function to adjust the colors of an L1 RGB (#76)

The atmosphere shifts things towards blue. The correction fixes that in
the composite to make it look nice.

doc/api/index.rst

@@ -26,6 +26,7 @@ Processing
     composite
     pansharpen
     equalize_histogram
+    adjust_l1_colors
     interpolate_missing
 
 Sample datasets

doc/plot-overlay.rst

@@ -40,7 +40,9 @@ have to show:
 .. jupyter-execute::
 
     # Make the composite
-    rgb = xls.composite(scene, rescale_to=(0, 0.2))
+    rgb = xls.composite(scene, rescale_to=(0, 0.15))
+    # Adjust the L1 colors to make it nicer and get rid of the blue glare
+    rgb = xls.adjust_l1_colors(rgb)
 
     # Plot the RGB and thermal separately
     fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))

xlandsat/__init__.py

@@ -3,7 +3,7 @@
 # SPDX-License-Identifier: MIT
 from . import datasets
 from ._composite import composite
-from ._enhancement import equalize_histogram
+from ._enhancement import adjust_l1_colors, equalize_histogram
 from ._interpolation import interpolate_missing
 from ._io import load_panchromatic, load_scene, save_scene
 from ._pansharpen import pansharpen

xlandsat/_enhancement.py

@@ -4,6 +4,7 @@
 """
 Operations to enhance composites. Basically wrappers for skimage.exposure.
 """
+import numpy as np
 import skimage.color
 import skimage.exposure
 
@@ -67,3 +68,48 @@ def equalize_histogram(composite, kernel_size=None, clip_limit=0.01):
         out_range=str(composite.values.dtype),
     )
     return result
+
+
+def adjust_l1_colors(composite, percentile=0.1):
+    """
+    Adjust the colors in an RGB composite from Level 1 data
+
+    Corrects the balance of the red, green, and blue bands in an RGB (true
+    color) composite made from Level 1 data (without atmospheric correction).
+    This is **not as accurate as atmospheric correction** but can lead to nicer
+    images in places where the atmospheric correction causes artifacts.
+
+    **Do not use the output for calculating indices.**
+
+    Parameters
+    ----------
+    composite : :class:`xarray.DataArray`
+        A composite, as generated by :func:`xlandsat.composite`.
+    percentile : float
+        The percentile range to use for the intensity rescaling. Will use the
+        ``percentile`` as the lower bound and ``100 - percentile`` for the
+        upper bound. Default is 0.1%.
+
+    Returns
+    -------
+    adjusted_composite : :class:`xarray.DataArray`
+        The composite after color adjustment, scaled back to the range of the
+        original composite.
+
+    Notes
+    -----
+
+    The correction raises each channel to the power of 1/2.2 (inverse of the
+    sRGB gamma function) and then rescales the transformed intensity to the
+    given percentile range.
+    """
+    output = composite.copy()
+    for i, channel in enumerate(["red", "green", "blue"]):
+        scaled = composite.sel(channel=channel).values ** (1 / 2.2)
+        vmin, vmax = np.percentile(scaled, (percentile, 100 - percentile))
+        output.values[:, :, i] = skimage.exposure.rescale_intensity(
+            scaled,
+            in_range=(vmin, vmax),
+            out_range=str(composite.values.dtype),
+        )
+    return output