Medoid Compositing

Seasonal Composite Landsat TM/ETM+ Images Using the Medoid (a Multi-Dimensional Median), Neil Flood, 2013, doi:10.3390/rs5126481



In [18]:

    
import ee
ee.Initialize()



In [19]:

    
from geetools import tools, composite, cloud_mask, indices



In [20]:

    
import ipygee as ui

Build a collection



In [21]:

    
p = ee.Geometry.Point(-72, -42)



In [22]:

    
col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\
        .filterBounds(p).filterDate('2017-01-01', '2017-12-01')\
        .map(cloud_mask.landsat8SRPixelQA())\
        .map(lambda img: img.addBands(indices.ndvi(img,'B5', 'B4')))\
        .limit(7)



In [23]:

    
ui.eprint(col.size())

Other simple composites to compare



In [24]:

    
max_ndvi = col.qualityMosaic('ndvi')



In [25]:

    
mosaic = col.mosaic()

Medoid



In [26]:

    
bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']

add date band before compositing



In [27]:

    
def add_date(img):
    date = tools.date.getDateBand(img)
    return img.addBands(date).copyProperties(date, ['day_since_epoch'])
col = col.map(add_date)



In [28]:

    
minval = ee.Number(col.aggregate_min('day_since_epoch'))



In [29]:

    
maxval = ee.Number(col.aggregate_max('day_since_epoch'))



In [30]:

    
medoid = composite.medoid(col, bands=bands)



In [31]:

    
medoid = medoid.set('min_date', minval, 'max_date', maxval)

Medoid without taking in count zero values



In [32]:

    
medoid_no_zeros = composite.medoid(col, bands=bands, discard_zeros=True)



In [33]:

    
medoid_no_zeros = medoid_no_zeros.set('min_date', minval, 'max_date', maxval)

Show on Map



In [35]:

    
Map = ui.Map()
Map.show()



In [36]:

    
vis = {'bands':['B5', 'B6','B4'], 'min':0, 'max':5000}



In [37]:

    
Map.addLayer(p)
Map.centerObject(p)



In [38]:

    
Map.addLayer(max_ndvi, vis, 'max NDVI')



In [39]:

    
Map.addLayer(mosaic, vis, 'simply Mosaic')



In [40]:

    
Map.addLayer(medoid, vis, 'Medoid')



In [41]:

    
Map.addLayer(medoid.select('date').randomVisualizer(), {'bands':['viz-red', 'viz-green', 'viz-blue'], 'min':0, 'max':255}, 'dates')



In [42]:

    
Map.addLayer(medoid_no_zeros, vis, 'Medoid without zero values')



In [43]:

    
Map.addLayer(medoid_no_zeros.select('date').randomVisualizer(), {'bands':['viz-red', 'viz-green', 'viz-blue'], 'min':0, 'max':255}, 'dates of medoid without zero values')

Extract data from images and compute locally to compare

Extract medoid values in point



In [44]:

    
medoid_values = tools.image.getValue(medoid.select(bands), p, scale=30, side='client')



In [45]:

    
medoid_values









    Out[45]:





{'B2': 77, 'B3': 188, 'B4': 105, 'B5': 2500, 'B6': 691, 'B7': 224}

List of values



In [46]:

    
medoid_values_list = [val for _, val in medoid_values.items()]



In [47]:

    
medoid_values_list









    Out[47]:





[77, 188, 105, 2500, 691, 224]

Extract values at point in each image of the collection



In [48]:

    
col_values = tools.imagecollection.getValues(col.select(bands), p, scale=30, side='client')

Get bandnames



In [49]:

    
col_key_list = []
for _, d in col_values.items():
    keys = []
    for k, v in d.items():
        keys.append(k)        
    col_key_list.append(keys)



In [50]:

    
col_key_list









    Out[50]:





[['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
 ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']]

Get values as a list



In [51]:

    
col_values_list = []
for _, d in col_values.items():
    values = []
    for _, v in d.items():
        if v:
            values.append(v)
        else:
            values.append(0)
    col_values_list.append(values)



In [52]:

    
col_values_list









    Out[52]:





[[0, 0, 0, 0, 0, 0],
 [107, 290, 159, 3142, 928, 307],
 [112, 272, 143, 3168, 870, 287],
 [87, 193, 107, 2465, 720, 245],
 [90, 210, 120, 2717, 813, 259],
 [77, 188, 105, 2500, 691, 224],
 [83, 200, 127, 2445, 698, 224]]

Medoid Method locally



In [53]:

    
def local_medoid(values):
    from copy import copy
    import math

    def distance(arr1, arr2):
        zipped = zip(arr1, arr2)
        accum = 0
        for a, b in zipped:
            calc = (a-b)*(a-b)
            accum += calc
        return math.sqrt(accum)

    def med(values):
        results = {}
        for i, val in enumerate(values):
            val = list(val)
            cop = copy(values)
            cop = [list(a) for a in cop]
            cop.remove(val)
            dist = 0
            for r in cop:
                r = list(r)
                d = distance(val, r)
                dist += d
            results[i] = dist

        return results
    
    def getmin(d):
        minval = min(d.values())
        for k, v in d.items():
            if v == minval:
                return k
    
    values = med(values)
    min_value = getmin(values)
    
    # return the index of the minimized sum as first argument, and all options as second
    return min_value, values

Compute medoid locally and compare



In [54]:

    
local = local_medoid(col_values_list)



In [55]:

    
local









    Out[55]:





(5,
 {0: 17251.127189606195,
  1: 5996.407209928899,
  2: 6022.895734074118,
  3: 4399.422196724633,
  4: 4593.00466857401,
  5: 4380.031017159888,
  6: 4461.227704323095})

Get the values that correspond to the medoid



In [56]:

    
min_values = col_values_list[local[0]]



In [57]:

    
min_values









    Out[57]:





[77, 188, 105, 2500, 691, 224]

Match bands with values



In [58]:

    
local_medoid = dict(zip(col_key_list[0], min_values))



In [59]:

    
local_medoid









    Out[59]:





{'B2': 77, 'B3': 188, 'B4': 105, 'B5': 2500, 'B6': 691, 'B7': 224}



In [60]:

    
medoid_values









    Out[60]:





{'B2': 77, 'B3': 188, 'B4': 105, 'B5': 2500, 'B6': 691, 'B7': 224}

Finally, compare values from medoid mosaic against locally computed medoid (from images values)



In [61]:

    
medoid_values == local_medoid









    Out[61]:





True



In [ ]: