Copyright © 2017 Naru. Tsutsumida (naru@kais.kyoto-u.ac.jp)
Using Google Earth Engine, Annually averaged MODIS MOD13Q1 NDVI data (250m spatial resolution, 16 days intervals) at a global scale are summarised during 2001-2016, then timeseries trend is calculated by Mann-kendall analysis.
Cloud masking and interpolation is not applied but code can be shown.
In [1]:
from IPython.display import Image, display, HTML
%matplotlib inline
from pylab import *
import datetime
import math
import time
import ee
ee.Initialize()
MODIS QA mask to filter out low quality pixels. several approaches for QA masks can be applied (but used only maskSummaryQA function)
In [2]:
def getQABits(image, start, end, newName):
#Compute the bits we need to extract.
p = 0
for i in range(start,(end+1)):
p += math.pow(2, i)
# Return a single band image of the extracted QA bits, giving the band
# a new name.
return image.select([0], [newName])\
.bitwiseAnd(p)\
.rightShift(start)
#A function to mask out cloudy pixels.
def maskClouds(img):
# Select the QA band.
QA = img.select('DetailedQA')
# Get the MOD_LAND_QA bits
internalCloud = getQABits(QA, 0, 1, 'MOD_LAND_QA')
# Return an image masking out cloudy areas.
return img.mask(internalCloud.eq(0))
##originally function for landsat
#https://groups.google.com/forum/#!searchin/google-earth-engine-developers/python$20bitwiseAnd%7Csort:relevance/google-earth-engine-developers/OYuUMjFr0Gg/GGtYWh4CAwAJ
def maskBadData(image):
invalid = image.select('DetailedQA').bitwiseAnd(0x6).neq(0)
clean = image.mask(invalid.Not())
return(clean)
def maskSummaryQA(img):
QA = img.select('SummaryQA').eq(0)
best = img.mask(QA)
return(best)
# function to add system time band
def addTimeBand(image):
return image.addBands(image.metadata('system:time_start').rename(["time"]))
in this demonstration the quality assurance (QA) mask is applied. Only the best quality pixel (QA=0) is picked up.
In [3]:
modisNDVI = ee.ImageCollection('MODIS/006/MOD13Q1') \
.select(['NDVI', "SummaryQA"]) \
.filterDate('2001-01-01', '2016-12-31') \
.sort('system:time_start')
count = modisNDVI.size().getInfo() ## total number of selected images
print('count:', count)
In [4]:
filteredMODIS = modisNDVI \
.map(maskSummaryQA) \
.select('NDVI')
toList(1,X). X=0,1,... first image as X=0, second as X=1....
check the info of first image
In [5]:
img1 = ee.Image(filteredMODIS.toList(1,0).get(0))
scale = img1.projection().nominalScale().getInfo()
props = img1.getInfo()['properties']
date = props['system:time_start']
system_time = datetime.datetime.fromtimestamp((date / 1000) - 3600)
date_str = system_time.strftime("%Y_%m_%d")
img1 = img1.set('bands', date_str)
##check metadata
print('scale:', scale) ##spatial resolution
print('DATE:', date_str) ##first date
In [6]:
## This field contains UNIX time in milliseconds.
timeField = 'system:time_start'
join = ee.Join.saveAll('images')
interval = 72 ##72 days
##ee.Filter.maxDifference:
##Creates a unary or binary filter that passes if the left and right operands, both numbers, are within a given maximum difference. If used as a join condition, this numeric difference is used as a join measure.
diffFilter = ee.Filter.maxDifference(difference = (1000 * 60 * 60 * 24) * interval,
leftField = timeField,
rightField = timeField)
NeighborJoin = join.apply(primary = filteredMODIS,
secondary = filteredMODIS,
condition = diffFilter)
def smooth_func(image):
collection = ee.ImageCollection.fromImages(image.get('images'))
return ee.Image(image).addBands(collection.mean().rename(['smooth']))
smoothed = ee.ImageCollection(NeighborJoin.map(smooth_func))
See a map of whole mean NDVI in 2001-2016.
Note that all data is in GEE server, not in this running environment.
In [8]:
regionfilter = ee.Geometry.Polygon([-170, 80, 0, 80, 170, 80, 170, -80, 10, -80, -170, -80]).toGeoJSON()
ndvi_palette = 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400, 3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
vizParams = {'min': -2000,
'max': 10000,
'region':regionfilter,
'palette': ndvi_palette}
img = smoothed.select('smooth').mean()
%config InlineBackend.figure_format = 'retina'
print(img.getThumbUrl(vizParams))
Image(url=img.getThumbUrl(vizParams), width=900, unconfined=True)
Out[8]:
After annual averaged NDVI datasets in 2001-2016 are calculated. Here Mann-Kendall trend test is applied.
Keep in mind as p-values are not calculated at this time on GEE (1st July 2017),
The statistically significant test is not available.
In [9]:
#create name list
yr = range(2001, 2017)
yr = map(str, yr)
num=len(range(2001, 2017))
yy = np.array(["Y"]*num)
years = np.core.defchararray.add(yy, yr)
st = np.array(["-01-01"]*num)
ed = np.array(["-12-31"]*num)
starts = np.core.defchararray.add(yr, st)
ends = np.core.defchararray.add(yr, ed)
In [10]:
y = 0
MODcoll = smoothed \
.filterDate(starts[y], ends[y]) \
.sort('system:time_start') \
.select('smooth')
start = starts[y]
end = ends[y]
avg = MODcoll.mean().rename([years[y]])
for y in range(1, 16):
MODcoll = smoothed \
.filterDate(starts[y], ends[y]) \
.sort('system:time_start') \
.select('smooth')
start = starts[y]
end = ends[y]
average = MODcoll.mean()
avg = avg.addBands(average.rename([years[y]]))
In [11]:
info_bands = avg.getInfo()['bands']
#print('Dimensions:', info_bands[0]['dimensions'])
print('Number of bands:', len(info_bands))
##see band names
for ids in range(0,len(info_bands),1):
print(info_bands[ids]['id'])
In [12]:
mk_ans = avg.reduce(ee.Reducer.kendallsCorrelation(1))
In [13]:
info_bands = mk_ans.getInfo()['bands']
#print('Dimensions:', info_bands[0]['dimensions'])
print('Number of bands:', len(info_bands))
##see bands
for ids in range(0,len(info_bands),1):
print(info_bands[ids]['id'])
In [14]:
RdBu_palette = '#B2182B, #D6604D, #F4A582, #FDDBC7, #F7F7F7, #D1E5F0, #92C5DE, #4393C3, #2166AC'
mk_tau = ee.Image(mk_ans.select('tau')).multiply(10000).int16()
url = mk_tau.getThumbUrl({
'min':-10000,
'max':10000,
'region':regionfilter,
'crs': 'EPSG:4326',
'palette': RdBu_palette
})
%config InlineBackend.figure_format = 'retina'
print(url)
Image(url=url, width=900, unconfined=True)
Out[14]:
Export the outputs to Google Drive.
It takes some time.
see: https://stackoverflow.com/questions/39219705/how-to-download-images-using-google-earth-engines-python-api
In [20]:
globalfilter = ee.Geometry.Polygon([-170, 80, 0, 80, 170, 80, 170, -80, 10, -80, -170, -80])
globalfilter = globalfilter['coordinates'][0]
task_config = {
'description': 'imageToDriveExample',
'scale': 231.65635826395828,
'region': globalfilter,
'maxPixels': 5e10
}
task = ee.batch.Export.image(mk_tau, 'kenn_0116', task_config)
task.start()
In [21]:
time.sleep(10)
while task.status()['state'] == 'RUNNING':
print 'Running...'
time.sleep(100)
print 'Done.', task.status()
In [ ]: