Revision ba7057c4
Added by Adam Wilson almost 11 years ago
| climate/procedures/ee_compile.R | ||
|---|---|---|
| 1 | 1 |
### Script to compile the monthly cloud data from earth engine into a netcdf file for further processing |
| 2 | 2 |
|
| 3 |
setwd("~/acrobates/adamw/projects/cloud")
|
|
| 4 | 3 |
library(raster) |
| 5 | 4 |
library(doMC) |
| 6 | 5 |
library(multicore) |
| 7 | 6 |
library(foreach) |
| 8 | 7 |
#library(doMPI) |
| 9 |
registerDoMC(10)
|
|
| 8 |
registerDoMC(4)
|
|
| 10 | 9 |
#beginCluster(4) |
| 11 | 10 |
|
| 12 |
tempdir="tmp" |
|
| 13 |
if(!file.exists(tempdir)) dir.create(tempdir) |
|
| 14 |
|
|
| 15 |
|
|
| 16 |
## Load list of tiles |
|
| 17 |
tiles=read.table("tile_lat_long_10d.txt",header=T)
|
|
| 18 |
|
|
| 19 |
### Build Tiles |
|
| 20 |
|
|
| 21 |
## bin sizes |
|
| 22 |
ybin=30 |
|
| 23 |
xbin=30 |
|
| 24 |
|
|
| 25 |
tiles=expand.grid(ulx=seq(-180,180-xbin,by=xbin),uly=seq(90,-90+ybin,by=-ybin)) |
|
| 26 |
tiles$h=factor(tiles$ulx,labels=paste("h",sprintf("%02d",1:length(unique(tiles$ulx))),sep=""))
|
|
| 27 |
tiles$v=factor(tiles$uly,labels=paste("v",sprintf("%02d",1:length(unique(tiles$uly))),sep=""))
|
|
| 28 |
tiles$tile=paste(tiles$h,tiles$v,sep="") |
|
| 29 |
tiles$urx=tiles$ulx+xbin |
|
| 30 |
tiles$ury=tiles$uly |
|
| 31 |
tiles$lrx=tiles$ulx+xbin |
|
| 32 |
tiles$lry=tiles$uly-ybin |
|
| 33 |
tiles$llx=tiles$ulx |
|
| 34 |
tiles$lly=tiles$uly-ybin |
|
| 35 |
tiles$cy=(tiles$uly+tiles$lry)/2 |
|
| 36 |
tiles$cx=(tiles$ulx+tiles$urx)/2 |
|
| 37 |
tiles=tiles[,c("tile","h","v","ulx","uly","urx","ury","lrx","lry","llx","lly","cx","cy")]
|
|
| 38 |
|
|
| 39 |
jobs=expand.grid(tile=tiles$tile,year=2000:2012,month=1:12) |
|
| 40 |
jobs[,c("ulx","uly","urx","ury","lrx","lry","llx","lly")]=tiles[match(jobs$tile,tiles$tile),c("ulx","uly","urx","ury","lrx","lry","llx","lly")]
|
|
| 41 |
|
|
| 42 |
## drop Janurary 2000 from list (pre-modis) |
|
| 43 |
jobs=jobs[!(jobs$year==2000&jobs$month==1),] |
|
| 44 |
|
|
| 45 |
#jobs=jobs[jobs$month==1,] |
|
| 46 |
|
|
| 47 |
|
|
| 48 |
#jobs=jobs[jobs$month==7,] |
|
| 49 |
## Run the python downloading script |
|
| 50 |
#system("~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin -159 20 -154.5 18.5 -year 2001 -month 6 -region test")
|
|
| 51 |
i=1 |
|
| 52 |
#testtiles=c("h02v07","h02v06","h02v08","h03v07","h03v06","h03v08")
|
|
| 53 |
#todo=which(jobs$tile%in%testtiles) |
|
| 54 |
#todo=todo[1:3] |
|
| 55 |
#todo=1 |
|
| 56 |
todo=1:nrow(jobs) |
|
| 57 |
|
|
| 58 |
checkcomplete=T |
|
| 59 |
if(checkcomplete&exists("df")){ #if desired (and "df" exists from below) drop complete date-tiles
|
|
| 60 |
todo=which(!paste(jobs$tile,jobs$year,jobs$month)%in%paste(df$region,df$year,df$month)) |
|
| 61 |
} |
|
| 62 |
|
|
| 63 |
writeLines(paste("Tiling options will produce",nrow(tiles),"tiles and ",nrow(jobs),"tile-months. Current todo list is ",length(todo)))
|
|
| 64 |
|
|
| 65 |
foreach(i=todo) %dopar%{
|
|
| 66 |
system(paste("python ~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin ",
|
|
| 67 |
jobs$ulx[i]," ",jobs$uly[i]," ",jobs$urx[i]," ",jobs$ury[i]," ",jobs$lrx[i]," ",jobs$lry[i]," ",jobs$llx[i]," ",jobs$lly[i]," ", |
|
| 68 |
" -year ",jobs$year[i]," -month ",jobs$month[i]," -region ",jobs$tile[i],sep="")) |
|
| 69 |
} |
|
| 11 |
wd="~/acrobates/adamw/projects/cloud" |
|
| 12 |
setwd(wd) |
|
| 70 | 13 |
|
| 71 | 14 |
|
| 72 | 15 |
## Get list of available files |
| ... | ... | |
| 87 | 30 |
table(df$year,df$month) |
| 88 | 31 |
|
| 89 | 32 |
## drop some if not complete |
| 90 |
df=df[df$month==1&df$year%in%c(2001:2002,2005),]
|
|
| 91 |
rerun=T # set to true to recalculate all dates even if file already exists
|
|
| 33 |
#df=df[df$month%in%1:9&df$year%in%c(2001:2012),]
|
|
| 34 |
rerun=F # set to true to recalculate all dates even if file already exists
|
|
| 92 | 35 |
|
| 93 | 36 |
## Loop over existing months to build composite netcdf files |
| 94 | 37 |
foreach(date=unique(df$date)) %dopar% {
|
| 95 |
## get date |
|
| 38 |
## get date
|
|
| 96 | 39 |
print(date) |
| 97 | 40 |
## Define output and check if it already exists |
| 98 |
tffile=paste(tempdir,"/mod09_",date,".tif",sep="")
|
|
| 41 |
vrtfile=paste(tempdir,"/mod09_",date,".vrt",sep="")
|
|
| 99 | 42 |
ncfile=paste(tempdir,"/mod09_",date,".nc",sep="") |
| 100 |
if(!rerun&file.exists(ncfile)) next |
|
| 43 |
tffile=paste(tempdir,"/mod09_",date,".tif",sep="") |
|
| 44 |
|
|
| 45 |
if(!rerun&file.exists(ncfile)) return(NA) |
|
| 101 | 46 |
## merge regions to a new netcdf file |
| 102 |
# system(paste("gdal_merge.py -tap -init 5 -o ",ncfile," -n -32768.000 -of netCDF -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))
|
|
| 103 |
system(paste("gdal_merge.py -o ",tffile," -init -32768 -n -32768.000 -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))
|
|
| 104 |
system(paste("gdal_translate -of netCDF ",tffile," ",ncfile," -ot Int16 "))
|
|
| 105 |
file.remove(tffile) |
|
| 47 |
# system(paste("gdal_merge.py -o ",tffile," -init -32768 -n -32768.000 -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))
|
|
| 48 |
system(paste("gdalbuildvrt -overwrite -srcnodata -32768 ",vrtfile," ",paste(df$path[df$date==date],collapse=" ")))
|
|
| 49 |
## Warp to WGS84 grid and convert to netcdf |
|
| 50 |
ops="-t_srs 'EPSG:4326' -multi -r cubic -te -90 -90 0 90 -tr 0.008333333333333 -0.008333333333333" |
|
| 51 |
ops="-t_srs 'EPSG:4326' -multi -r cubic -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333" |
|
| 52 |
|
|
| 53 |
system(paste("gdalwarp -overwrite ",ops," -srcnodata -32768 -dstnodata -32768 -of netCDF ",vrtfile," ",ncfile," -ot Int16"))
|
|
| 54 |
# system(paste("gdalwarp -overwrite ",ops," -srcnodata -32768 -dstnodata -32768 -of netCDF ",vrtfile," ",tffile," -ot Int16"))
|
|
| 55 |
|
|
| 56 |
setwd(wd) |
|
| 106 | 57 |
system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
| 107 | 58 |
## create temporary nc file with time information to append to MOD06 data |
| 108 | 59 |
cat(paste("
|
| ... | ... | |
| 122 | 73 |
## add other attributes |
| 123 | 74 |
system(paste("ncrename -v Band1,CF ",ncfile,sep=""))
|
| 124 | 75 |
system(paste("ncatted ",
|
| 125 |
" -a units,CF,o,c,\"Proportion Days Cloudy\" ",
|
|
| 76 |
" -a units,CF,o,c,\"%\" ",
|
|
| 126 | 77 |
# " -a valid_range,CF,o,b,\"0,100\" ", |
| 127 | 78 |
" -a scale_factor,CF,o,f,\"0.1\" ", |
| 128 |
" -a long_name,CF,o,c,\"Proportion cloudy days (%)\" ", |
|
| 79 |
" -a _FillValue,CF,o,f,\"-32768\" ", |
|
| 80 |
" -a long_name,CF,o,c,\"Cloud Frequency(%)\" ", |
|
| 81 |
" -a NETCDF_VARNAME,CF,o,c,\"Cloud Frequency(%)\" ", |
|
| 82 |
" -a title,global,o,c,\"Cloud Climatology from MOD09 Cloud Mask\" ", |
|
| 83 |
" -a institution,global,o,c,\"Jetz Lab, EEB, Yale University, New Haven, CT\" ", |
|
| 84 |
" -a source,global,o,c,\"Derived from MOD09GA Daily Data\" ", |
|
| 85 |
" -a comment,global,o,c,\"Developed by Adam M. Wilson (adam.wilson@yale.edu / http://adamwilson.us)\" ", |
|
| 129 | 86 |
ncfile,sep="")) |
| 130 | 87 |
|
| 131 | 88 |
## add the fillvalue attribute back (without changing the actual values) |
| ... | ... | |
| 137 | 94 |
} |
| 138 | 95 |
print(paste(basename(ncfile)," Finished")) |
| 139 | 96 |
|
| 140 |
} |
|
| 141 |
|
|
| 142 |
|
|
| 143 | 97 |
|
| 98 |
} |
|
| 144 | 99 |
|
| 145 | 100 |
|
| 146 | 101 |
### merge all the tiles to a single global composite |
| 147 | 102 |
#system(paste("ncdump -h ",list.files(tempdir,pattern="mod09.*.nc$",full=T)[10]))
|
| 148 |
system(paste("cdo -O mergetime ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/mod09.nc"))
|
|
| 103 |
system(paste("cdo -O mergetime ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_daily.nc"))
|
|
| 149 | 104 |
|
| 105 |
# Overall mean |
|
| 106 |
system(paste("cdo -O -chname,CF,CF_annual -timmean data/cloud_daily.nc data/cloud_mean.nc"))
|
|
| 150 | 107 |
|
| 151 | 108 |
### generate the monthly mean and sd |
| 152 | 109 |
#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))
|
| 153 |
xosystem(paste("cdo -O -ymonmean data/mod09.nc data/mod09_clim_mean.nc"))
|
|
| 154 |
system(paste("cdo -O -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim_sd.nc"))
|
|
| 110 |
system(paste("cdo -O -ymonmean data/cloud_daily.nc data/cloud_ymonmean.nc"))
|
|
| 111 |
system(paste("cdo -O -chname,CF,CF_sd -ymonstd data/cloud_daily.nc data/cloud_ymonsd.nc"))
|
|
| 112 |
|
|
| 113 |
#if(!file.exists("data/mod09_metrics.nc")) {
|
|
| 114 |
# system("cdo -chname,CF,CFmin -timmin data/mod09_clim_mean.nc data/mod09_min.nc")
|
|
| 115 |
# system("cdo -chname,CF,CFmax -timmax data/mod09_clim_mean.nc data/mod09_max.nc")
|
|
| 116 |
# system("cdo -chname,CF,CFsd -timstd data/mod09_clim_mean.nc data/mod09_std.nc")
|
|
| 117 |
# system("cdo -f nc2 merge data/mod09_std.nc data/mod09_min.nc data/mod09_max.nc data/mod09_metrics.nc")
|
|
| 118 |
system("cdo merge -chname,CF,CFmin -timmin data/cloud_clim_mean.nc -chname,CF,CFmax -timmax data/cloud_clim_mean.nc -chname,CF,CFsd -timstd data/cloud_clim_mean.nc data/cloud_metrics.nc")
|
|
| 119 |
#} |
|
| 120 |
|
|
| 121 |
|
|
| 122 |
|
|
| 123 |
|
|
| 124 |
|
|
| 155 | 125 |
|
| 156 |
# Overall mean |
|
| 157 |
system(paste("cdo -O -chname,CF,CF_annual -timmean data/mod09.nc data/mod09_clim_mac.nc"))
|
|
| 158 | 126 |
|
| 159 | 127 |
### Long term summaries |
| 160 | 128 |
seasconc <- function(x,return.Pc=T,return.thetat=F) {
|
| ... | ... | |
| 187 | 155 |
} |
| 188 | 156 |
|
| 189 | 157 |
|
| 158 |
|
|
| 190 | 159 |
## read in monthly dataset |
| 191 | 160 |
mod09=brick("data/mod09_clim_mean.nc",varname="CF")
|
| 192 | 161 |
plot(mod09[1]) |
Also available in: Unified diff
Separated ee download from compile script. Switched to use of vrt as intermediate file form