diff --git a/doc/api/index.rst b/doc/api/index.rst index c89f923..5a69dab 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -26,6 +26,7 @@ Processing composite pansharpen equalize_histogram + adjust_l1_colors interpolate_missing Sample datasets diff --git a/doc/plot-overlay.rst b/doc/plot-overlay.rst index f16b91c..4bdbf4e 100644 --- a/doc/plot-overlay.rst +++ b/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)) diff --git a/xlandsat/__init__.py b/xlandsat/__init__.py index 93d4d90..0748713 100644 --- a/xlandsat/__init__.py +++ b/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 diff --git a/xlandsat/_enhancement.py b/xlandsat/_enhancement.py index b0153ff..1188f3e 100644 --- a/xlandsat/_enhancement.py +++ b/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