Histogram Matching

Tutorial created by **Aaron Zuspan**: GitHub | Twitter

GitHub Repo: https://github.com/davemlz/eemont
PyPI link: https://pypi.org/project/eemont/
Conda-forge: https://anaconda.org/conda-forge/eemont
Documentation: https://eemont.readthedocs.io/
More tutorials: https://github.com/davemlz/eemont/tree/master/docs/tutorials

Let’s start!

If required, please uncomment:

[1]:

# !pip install eemont
# !pip install geemap

Import the required packages.

[2]:

import ee, eemont, geemap

Authenticate and Initialize Earth Engine and geemap.

[3]:

Map = geemap.Map()

Matching Landsat 5 to Landsat 8

Histogram matching performs band-wise adjustments to match the spectral characteristics of one image to an overlapping image. It can be used for images taken at different times or by different sensors to allow comparisons.

In this example we’ll match a Landsat 5 image to a Landsat 8 image. We’ll use images from different seasons to demonstrate a dramatic example of histogram adjustment, but note that significant changes between images such as land cover, seasonality, or clouds and snow may create inaccurate imagery.

Load and visualize a Landsat 5 image from winter of 2011.

[4]:

source = ee.Image("LANDSAT/LT05/C01/T1_SR/LT05_195028_20110208")

Map.addLayer(source, {"min": 170, "max": 2110, "bands": ["B3", "B2", "B1"]}, "Landsat 5 - Source")
pt = ee.Geometry.Point([6.6346, 46.6989])
Map.centerObject(pt, zoom=12)
Map

Load a Landsat 8 image from the same location in summer of 2018. Notice how the Landsat 8 image is much greener.

[5]:

target = ee.Image("LANDSAT/LC08/C02/T1_L2/LC08_196027_20130714")
Map.addLayer(target, {"min": 7800, "max": 13700, "bands": ["SR_B4", "SR_B3", "SR_B2"]}, "Landsat 8 - Target")
Map

Specify a region to match histograms within. Image histograms will be calculated within this region, so choose an area that is representative of each image and that does not contain clouds. If no geometry is provided, the full source image will be used. We’ll visualize the region later using an outline.

[6]:

region = ee.Geometry.Polygon(
    [[[6.499245, 46.640447],
      [6.499245, 46.763351],
      [6.779451, 46.763351],
      [6.779451, 46.640447],
      [6.499245, 46.640447]]]
)

regionOutline = ee.Image().byte().paint(featureCollection=region, color=1, width=1)

Specify the corresponding band pairs to match using a dictionary with the source bands as keys and the target bands as values. Only the bands specified will be matched and included in the output image. If you’re not sure which bands correspond between platforms, try checking the band descriptions in the Earth Engine Data Catalog or use `ee.Image.getSTAC() <https://colab.research.google.com/github/davemlz/eemont/blob/master/tutorials/019-Checking-STAC-Info.ipynb>`__.

[7]:

bands = {
    # Red
    "B3": "SR_B4",
    # Green
    "B2": "SR_B3",
    # Blue
    "B1": "SR_B2"
}

Histogram-match the Landsat 5 image to the Landsat 8 image and visualize the result. Notice how the adjusted Landsat 5 image is much greener, matching the Landsat 8 image.

[8]:

matched = source.matchHistogram(target, bands, geometry=region)
Map.addLayer(matched, {"min": 7800, "max": 13700, "bands": ["B3", "B2", "B1"]}, "Landsat 5 - Matched")
Map.addLayer(regionOutline, {"palette": "#00FFFF"}, "Region")
Map

Matching Sentinel-2 to MODIS

Create a fresh map.

[9]:

Map = geemap.Map()

In this example, we’ll match the histogram of a Sentinel-2 image to a MODIS image taken on a similar date to minimize atmospheric effects in the Sentinel image. First, we’ll load the MODIS imagery.

[10]:

target = ee.Image("MODIS/006/MOD09A1/2018_08_05")
Map.addLayer(target, {"min": 600, "max": 4000, "bands": ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"]}, "MODIS - Source")
pt = ee.Geometry.Point([29.4450, 12.9945])
Map.centerObject(pt, zoom=9)
Map

Next, we’ll load a Sentinel-2 image. Notice how the colors do not match the MODIS image.

[11]:

source = ee.Image("COPERNICUS/S2/20180923T081641_20180923T083023_T35PQQ")
Map.addLayer(source, {"min": 600, "max": 4000, "bands": ["B4", "B3", "B2"]}, "Sentinel-2 - Source")
Map

Specify which band pairs to match.

[12]:

bands = {
    # Red
    "B4": "sur_refl_b01",
    # Green
    "B3": "sur_refl_b04",
    # Blue
    "B2": "sur_refl_b03"
}

Histogram-match the Sentinel-2 image to the MODIS image. Because we’re not specifying a region, the entire Sentinel-2 image footprint will be used for matching.

Notice how the matched image now blends in with the MODIS image.

[13]:

matched = source.matchHistogram(target, bands)
Map.addLayer(matched, {"min": 600, "max": 4000, "bands": ["B4", "B3", "B2"]}, "Sentinel-2 - Matched")
Map

Time Series By Regions (with Conversion to Pandas and Visualization with Seaborn)

Panchromatic Sharpening