Skip to content

Instantly share code, notes, and snippets.

@flacle
Last active November 1, 2021 15:39
Show Gist options
  • Save flacle/fac881e4d9348667b2729653be940c4f to your computer and use it in GitHub Desktop.
Save flacle/fac881e4d9348667b2729653be940c4f to your computer and use it in GitHub Desktop.
Histogram match - Google Earth Engine
/**
* # Histogram matching
* # Google Earth Engine with S1 SAR
* # A very small experiment :D
* # Author: Francis Laclé
* # Year: 2019
*/
/**
* Function that returns a normalized image using unitScale()
* @param image the image object
* @param band the band to select
* @param aoi the area of interest for scoping the reducer
*/
function getNormalizedImage(image, band, aoi) {
var redux = getReducedImage(image, ee.Reducer.minMax(), aoi);
return image.unitScale(
ee.Number(redux.get(band + '_min')),
ee.Number(redux.get(band + '_max'))
);
}
/**
* Function wrapper for a reducer
* @param image the image object
* @param reducer the reducer to apply, see API docs
* @param aoi the area of interest for scoping the reducer
*/
function getReducedImage(image, reducer, aoi) {
return image.reduceRegion({
reducer: reducer,
geometry: aoi,
scale: 10,
bestEffort: true
});
}
/**
* Function to print a histogram in the range of 0..255
* @param y is the y-value for the chart, in this case the frequencies
* @param label the label used for the series
*/
function printHistogram(y, label) {
var x = ee.List.sequence(0, 255);
var chart = ui.Chart.array.values(y, 0, x)
.setSeriesNames(label)
.setOptions({
title: 'Histogram',
hAxis: {'title': 'Bin'},
vAxis: {'title': 'Frequency'},
pointSize: 3,
series: {
0: {color: 'ff0000'},
1: {color: '00ff00'},
2: {color: '0000ff'}
}
});
print(chart);
}
/**
* Function that returns a histogram array object
* @param image the image object
* @param band the band to select
* @param aoi the area of interest for scoping the reducer
*/
function getHistogram(image, band, aoi) {
var redux = image.reduceRegion({
reducer: ee.Reducer.histogram(256, 1),
geometry: aoi,
scale: 10,
maxPixels: 1e9
});
var histoDict = ee.Dictionary(redux.get(band)).get('histogram');
return histoDict;
}
/**
* Function that returns a cumulative distribution function
* as a list object
* @param image the image object
* @param band the band to select
* @param aoi the area of interest for scoping the reducer
*/
function getCDF(image, band, aoi) {
var histoArr = ee.List(getHistogram(image.multiply(255), band, aoi));
var len = histoArr.length().subtract(1);
var nPixels = ee.Array(histoArr).accum(0).toList().get(len);
var CDF = ee.Array(histoArr).accum(0).divide(nPixels).multiply(len);
return [CDF.round().toList(), len];
}
/**
* Function that returns a source image that has been
* matched with the target's image histogram using
* the remap function from Earth Engine
* @param imageR the source image object
* @param imageZ the target image object
* @param bandR the source band to select
* @param bandZ the target band to select
* @param aoi the area of interest for scoping the reducer
*/
function matchHistogram(imageR, imageZ, bandR, bandZ, aoi) {
var cdfR = getCDF(imageR, bandR, aoi);
var cdfZ = getCDF(imageZ, bandZ, aoi);
var lenR = cdfR[1];
var lenZ = cdfZ[1];
var map = [];
var cdfR0 = cdfR[0].getInfo();
var cdfZ0 = cdfZ[0].getInfo();
for (var r = 0; r < cdfR0.length; r++) {
var max = Infinity;
var min = 0;
for (var z = 0; z < cdfZ0.length; z++) {
var diff = Math.pow(cdfZ0[z] - cdfR0[r], 2);
if (diff < max) {
max = diff;
min = z;
}
}
map.push(min);
}
var seq = ee.List.sequence(0, lenR);
var imageM = imageR.multiply(lenR).round().remap(seq, map, 0).divide(lenR);
return imageM.select(['remapped'], [bandR]);
}
// Aruba ~midpoint: -69.96542164802798,12.51783066594331
var poi = ee.Geometry.Point(-69.96542164802798,12.51783066594331);
var bounds = 12000;
var aoi = poi.buffer(bounds).bounds();
var zoom = 12;
Map.centerObject(aoi, zoom);
var img1 = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterDate('2016-01-01', '2016-06-01')
.sort('system:time_start', true)
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.eq('resolution', 'H'))
.filter(ee.Filter.eq('resolution_meters', 10))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.select('VV')
.filterBounds(aoi);
var img2 = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterDate('2016-06-01', '2016-12-01')
.sort('system:time_start', true)
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.eq('resolution', 'H'))
.filter(ee.Filter.eq('resolution_meters', 10))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.select('VV')
.filterBounds(aoi);
img1 = getNormalizedImage(img1.mosaic(), 'VV', aoi);
img2 = getNormalizedImage(img2.mosaic(), 'VV', aoi);
var match = matchHistogram(img1, img2, 'VV', 'VV', aoi);
Map.addLayer(img1.clip(aoi), {}, 'img1');
Map.addLayer(img2.clip(aoi), {}, 'img2');
Map.addLayer(match.clip(aoi), {}, 'match');
var vis = {
min: 0.0,
max: 1.0,
bands: ['VV'],
};
print(img1.clip(aoi).getThumbURL(vis));
print(img2.clip(aoi).getThumbURL(vis));
print(match.clip(aoi).getThumbURL(vis));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment