1
|
## Example script that downloads data from Google Earth Engine using the python API
|
2
|
## MODIS MOD09GA data is processed to extract the MOD09 cloud flag and calculate monthly cloud frequency
|
3
|
|
4
|
|
5
|
## import some libraries
|
6
|
import ee
|
7
|
from ee import mapclient
|
8
|
import ee.mapclient
|
9
|
import datetime
|
10
|
import wget
|
11
|
import os
|
12
|
from subprocess import call
|
13
|
|
14
|
#import logging
|
15
|
#logging.basicConfig()
|
16
|
|
17
|
## set working directory (where files will be downloaded)
|
18
|
os.chdir('/home/adamw/acrobates/adamw/projects/cloud/data/mod09')
|
19
|
|
20
|
MY_SERVICE_ACCOUNT = '511722844190@developer.gserviceaccount.com' # replace with your service account
|
21
|
MY_PRIVATE_KEY_FILE = '/home/adamw/EarthEngine-privatekey.p12' # replace with you private key file path
|
22
|
|
23
|
#MY_SERVICE_ACCOUNT = '205878743334-4mrtqgu0n5rnsv1vanrvv6atqk6vu8am@developer.gserviceaccount.com'
|
24
|
#MY_PRIVATE_KEY_FILE = '/home/adamw/EarthEngine_Jeremy-privatekey.p12'
|
25
|
|
26
|
ee.Initialize(ee.ServiceAccountCredentials(MY_SERVICE_ACCOUNT, MY_PRIVATE_KEY_FILE))
|
27
|
|
28
|
## set map center to speed up viewing
|
29
|
#ee.mapclient.centerMap(-121.767, 46.852, 11)
|
30
|
|
31
|
#///////////////////////////////////
|
32
|
#// Function to extract cloud flags
|
33
|
def getmod09(img): return(img.select(['state_1km']).expression("((b(0)/1024)%2)>0.5"));
|
34
|
# added the >0.5 because some values are coming out >1. Need to look into this further as they should be bounded 0-1...
|
35
|
|
36
|
#// Date ranges
|
37
|
yearstart=2000
|
38
|
yearstop=2012
|
39
|
monthstart=1
|
40
|
monthstop=12
|
41
|
|
42
|
#////////////////////////////////////////////////////
|
43
|
# Loop through months and get monthly % missing data
|
44
|
|
45
|
## set a year-month if you don't want to run the loop (for testing)
|
46
|
#year=2001
|
47
|
#month=2
|
48
|
|
49
|
|
50
|
## define the regions to be processed
|
51
|
regions=['[[-180, -60], [-180, 0], [0, 0], [0, -60]]', # SW
|
52
|
'[[-180, 0], [-180, 90], [0, 90], [0, 0]]', # NW
|
53
|
'[[0, 0], [0, 90], [180, 90], [180, 0]]', # NE
|
54
|
'[[0, 0], [0, -60], [180, -60], [180, 0]]'] # SE
|
55
|
## name the regions (these names will be used in the file names
|
56
|
## must be same length as regions list above
|
57
|
rnames=['SW','NW','NE','SE']
|
58
|
|
59
|
## Loop over regions, years, months to generate monthly timeseries
|
60
|
for r in range(0,len(regions)): # loop over regions
|
61
|
for year in range(yearstart,yearstop+1): # loop over years
|
62
|
for month in range(monthstart,monthstop+1): # loop over months
|
63
|
print('Processing '+rnames[r]+"_"+str(year)+'_'+str(month))
|
64
|
|
65
|
# output filename
|
66
|
filename='mod09_'+rnames[r]+"_"+str(year)+"_"+str(month)
|
67
|
|
68
|
# Check if file already exists and continue if so...
|
69
|
if(os.path.exists(filename+".zip")):
|
70
|
print("File exists:"+filename)
|
71
|
continue
|
72
|
|
73
|
# MOD09 internal cloud flag for this year-month
|
74
|
# to filter by a date range: filterDate(datetime.datetime(yearstart,monthstart,1),datetime.datetime(yearstop,monthstop,31))
|
75
|
mod09 = ee.ImageCollection("MOD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);
|
76
|
# calculate mean cloudiness (%), rename band, multiply by 100, and convert to integer
|
77
|
mod09a=mod09.mean().select([0], ['MOD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).byte();
|
78
|
|
79
|
# Next few lines for testing only
|
80
|
# print info to confirm there is data
|
81
|
# mod09a.getInfo()
|
82
|
# add to plot to confirm it's working
|
83
|
# ee.mapclient.addToMap(mod09a, {'range': '0,100'}, 'MOD09')
|
84
|
|
85
|
# build the URL and name the object (so that when it's unzipped we know what it is!)
|
86
|
path =mod09a.getDownloadUrl({
|
87
|
'name': filename, # name the file (otherwise it will be a uninterpretable hash)
|
88
|
'scale': 1000, # resolution in meters
|
89
|
'crs': 'EPSG:4326', # MODIS sinusoidal
|
90
|
'region': regions[r] # region defined above
|
91
|
});
|
92
|
|
93
|
# Sometimes EE will serve a corrupt zipped file with no error
|
94
|
# to check this, use a while loop that keeps going till there is an unzippable file.
|
95
|
# This has the potential for an infinite loop...
|
96
|
|
97
|
while(not(os.path.exists(filename+".zip"))):
|
98
|
# download with wget
|
99
|
print("Downloading "+filename)
|
100
|
wget.download(path)
|
101
|
# try to unzip it
|
102
|
print("Unzipping "+filename)
|
103
|
zipstatus=call("unzip "+filename+".zip",shell=True)
|
104
|
# if file doesn't exists or it didn't unzip, remove it and try again
|
105
|
if(zipstatus==9):
|
106
|
print("ERROR: "+filename+" unzip-able")
|
107
|
os.remove(filename)
|
108
|
|
109
|
print 'Finished'
|
110
|
|
111
|
|