Revision f5ef0596
Added by Adam Wilson almost 11 years ago
| climate/procedures/ee.MOD09.py | ||
|---|---|---|
| 1 |
#!/usr/bin/env python |
|
| 2 |
|
|
| 1 | 3 |
## Example script that downloads data from Google Earth Engine using the python API |
| 2 | 4 |
## MODIS MOD09GA data is processed to extract the MOD09 cloud flag and calculate monthly cloud frequency |
| 3 | 5 |
|
| ... | ... | |
| 9 | 11 |
import datetime |
| 10 | 12 |
import wget |
| 11 | 13 |
import os |
| 14 |
import sys |
|
| 12 | 15 |
from subprocess import call |
| 13 | 16 |
|
| 14 |
#import logging |
|
| 15 |
#logging.basicConfig() |
|
| 17 |
import logging |
|
| 18 |
logging.basicConfig(filename='error.log',level=logging.DEBUG) |
|
| 19 |
|
|
| 20 |
def Usage(): |
|
| 21 |
print('Usage: ee.MOD9.py -projwin ulx uly lrx lry -year year -month month -regionname 1')
|
|
| 22 |
sys.exit( 1 ) |
|
| 23 |
|
|
| 24 |
ulx = float(sys.argv[2]) |
|
| 25 |
uly = float(sys.argv[3]) |
|
| 26 |
lrx = float(sys.argv[4]) |
|
| 27 |
lry = float(sys.argv[5]) |
|
| 28 |
year = int(sys.argv[7]) |
|
| 29 |
month = int(sys.argv[9]) |
|
| 30 |
regionname = str(sys.argv[11]) |
|
| 31 |
|
|
| 32 |
#``` |
|
| 33 |
#ulx=-159 |
|
| 34 |
#uly=20 |
|
| 35 |
#lrx=-154.5 |
|
| 36 |
#lry=18.5 |
|
| 37 |
#year=2001 |
|
| 38 |
#month=6 |
|
| 39 |
#``` |
|
| 40 |
## Define output filename |
|
| 41 |
output=regionname+'_'+str(year)+'_'+str(month) |
|
| 16 | 42 |
|
| 17 | 43 |
## set working directory (where files will be downloaded) |
| 18 | 44 |
os.chdir('/mnt/data2/projects/cloud/mod09')
|
| ... | ... | |
| 20 | 46 |
MY_SERVICE_ACCOUNT = '511722844190@developer.gserviceaccount.com' # replace with your service account |
| 21 | 47 |
MY_PRIVATE_KEY_FILE = '/home/adamw/EarthEngine-privatekey.p12' # replace with you private key file path |
| 22 | 48 |
|
| 23 |
#MY_SERVICE_ACCOUNT = '205878743334-4mrtqgu0n5rnsv1vanrvv6atqk6vu8am@developer.gserviceaccount.com' |
|
| 24 |
#MY_PRIVATE_KEY_FILE = '/home/adamw/EarthEngine_Jeremy-privatekey.p12' |
|
| 25 |
|
|
| 26 | 49 |
ee.Initialize(ee.ServiceAccountCredentials(MY_SERVICE_ACCOUNT, MY_PRIVATE_KEY_FILE)) |
| 27 | 50 |
|
| 28 | 51 |
## set map center to speed up viewing |
| ... | ... | |
| 33 | 56 |
def getmod09(img): return(img.select(['state_1km']).expression("((b(0)/1024)%2)>0.5"));
|
| 34 | 57 |
# 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 | 58 |
|
| 36 |
#// Date ranges |
|
| 37 |
yearstart=2000 |
|
| 38 |
yearstop=2012 |
|
| 39 |
monthstart=1 |
|
| 40 |
monthstop=12 |
|
| 41 |
|
|
| 42 | 59 |
#//////////////////////////////////////////////////// |
| 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 |
#r=1 |
|
| 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 | 60 |
|
| 65 | 61 |
# output filename |
| 66 |
filename='mod09_'+rnames[r]+"_"+str(year)+"_"+str(month) |
|
| 67 |
unzippedfilename='mod09_'+rnames[r]+"_"+str(year)+"_"+str(month)+".MOD09_"+str(year)+"_"+str(month)+".tif" |
|
| 62 |
unzippedfilename=output+".mod09.tif" |
|
| 68 | 63 |
|
| 69 | 64 |
# Check if file already exists and continue if so... |
| 70 |
if(os.path.exists(unzippedfilename)): |
|
| 71 |
print("File exists:"+filename)
|
|
| 72 |
continue |
|
| 65 |
if(os.path.exists(unzippedfilename)): |
|
| 66 |
sys.exit("File exists:"+output)
|
|
| 73 | 67 |
|
| 74 |
# MOD09 internal cloud flag for this year-month |
|
| 68 |
|
|
| 69 |
##################################################### |
|
| 70 |
# Processing Function |
|
| 71 |
# MOD09 internal cloud flag for this year-month |
|
| 75 | 72 |
# to filter by a date range: filterDate(datetime.datetime(yearstart,monthstart,1),datetime.datetime(yearstop,monthstop,31)) |
| 76 |
mod09 = ee.ImageCollection("MOD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);
|
|
| 73 |
mod09 = ee.ImageCollection("MOD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);
|
|
| 77 | 74 |
# myd09 = ee.ImageCollection("MYD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);
|
| 78 | 75 |
# calculate mean cloudiness (%), rename band, multiply by 100, and convert to integer |
| 79 |
mod09a=mod09.mean().select([0], ['MOD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).byte(); |
|
| 80 |
# myd09a=myd09.mean().select([0], ['MYD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).byte(); |
|
| 81 |
|
|
| 82 |
``` |
|
| 76 |
mod09a=mod09.mean().select([0], ['mod09']).multiply(ee.Image(1000)).int16(); |
|
| 77 |
# myd09a=myd09.mean().select([0], ['MYD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).int8(); |
|
| 78 |
## Set data equal to whatver you want downloaded |
|
| 79 |
data=mod09a |
|
| 80 |
###################################################### |
|
| 81 |
|
|
| 82 |
## define region for download |
|
| 83 |
region=[ulx,lry], [ulx, uly], [lrx, uly], [lrx, lry] #h11v08 |
|
| 84 |
strregion=str(list(region)) |
|
| 83 | 85 |
# Next few lines for testing only |
| 84 | 86 |
# print info to confirm there is data |
| 85 |
mod09a.getInfo() |
|
| 86 |
myd09a.getInfo() |
|
| 87 |
#data.getInfo() |
|
| 88 |
|
|
| 89 |
## print a status update |
|
| 90 |
print(output+' Processing.... Coords:'+strregion) |
|
| 91 |
|
|
| 87 | 92 |
|
| 88 | 93 |
# add to plot to confirm it's working |
| 89 |
ee.mapclient.addToMap(mod09a, {'range': '0,100'}, 'MOD09')
|
|
| 90 |
ee.mapclient.addToMap(myd09a, {'range': '0,100'}, 'MOD09')
|
|
| 91 |
``` |
|
| 94 |
#ee.mapclient.addToMap(data, {'range': '0,100'}, 'MOD09')
|
|
| 95 |
#``` |
|
| 96 |
|
|
| 97 |
# TODO: |
|
| 98 |
# use MODIS projection |
|
| 92 | 99 |
|
| 93 | 100 |
# build the URL and name the object (so that when it's unzipped we know what it is!) |
| 94 |
path =mod09a.getDownloadUrl({
|
|
| 95 |
'name': filename, # name the file (otherwise it will be a uninterpretable hash)
|
|
| 96 |
'scale': 1000, # resolution in meters
|
|
| 97 |
'crs': 'EPSG:4326', # MODIS sinusoidal
|
|
| 98 |
'region': regions[r] # region defined above
|
|
| 99 |
});
|
|
| 101 |
path =mod09a.getDownloadUrl({
|
|
| 102 |
'name': output, # name the file (otherwise it will be a uninterpretable hash)
|
|
| 103 |
'scale': 926, # resolution in meters
|
|
| 104 |
'crs': 'EPSG:4326', # projection
|
|
| 105 |
'region': strregion # region defined above
|
|
| 106 |
}); |
|
| 100 | 107 |
|
| 101 | 108 |
# Sometimes EE will serve a corrupt zipped file with no error |
| 102 | 109 |
# to check this, use a while loop that keeps going till there is an unzippable file. |
| 103 | 110 |
# This has the potential for an infinite loop... |
| 104 | 111 |
|
| 105 |
while(not(os.path.exists(filename+".zip"))): |
|
| 106 |
# download with wget |
|
| 107 |
print("Downloading "+filename)
|
|
| 108 |
wget.download(path) |
|
| 112 |
#if(not(os.path.exists(output+".tif"))): |
|
| 113 |
# download with wget |
|
| 114 |
print("Downloading "+output)
|
|
| 115 |
wget.download(path) |
|
| 116 |
#call(["wget"+path,shell=T]) |
|
| 109 | 117 |
# try to unzip it |
| 110 |
print("Unzipping "+filename)
|
|
| 111 |
zipstatus=call("unzip "+filename+".zip",shell=True)
|
|
| 118 |
print("Unzipping "+output)
|
|
| 119 |
zipstatus=call("unzip "+output+".zip",shell=True)
|
|
| 112 | 120 |
# if file doesn't exists or it didn't unzip, remove it and try again |
| 113 |
if(zipstatus==9): |
|
| 114 |
print("ERROR: "+filename+" unzip-able")
|
|
| 115 |
os.remove(filename) |
|
| 116 |
|
|
| 117 |
print 'Finished' |
|
| 121 |
if(zipstatus==9): |
|
| 122 |
sys.exit("File exists:"+output)
|
|
| 123 |
# print("ERROR: "+output+" unzip-able")
|
|
| 124 |
# os.remove(output+".zip") |
|
| 125 |
|
|
| 126 |
## delete the zipped file (the unzipped version is kept) |
|
| 127 |
os.remove(output+".zip") |
|
| 128 |
|
|
| 129 |
print(output+' Finished!') |
|
| 118 | 130 |
|
| 119 | 131 |
|
Also available in: Unified diff
Updated EE script to run in parallel and minor edits to post-processing