Revision aba23d60
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_Climatology.r | ||
|---|---|---|
| 1 | 1 |
### Process a folder of daily MOD35 HDF files to produce a climatology |
| 2 | 2 |
|
| 3 |
.libPaths("/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
|
| 4 |
|
|
| 3 | 5 |
## import commandline arguments |
| 4 |
library(getopt) |
|
| 6 |
library(getopt,lib="/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
|
| 5 | 7 |
## get options |
| 6 | 8 |
opta <- getopt(matrix(c( |
| 7 | 9 |
'tile', 't', 1, 'character', |
| ... | ... | |
| 11 | 13 |
tile=opta$tile #tile="h11v08" |
| 12 | 14 |
verbose=opta$verbose #print out extensive information for debugging? |
| 13 | 15 |
|
| 16 |
## set working directory |
|
| 17 |
setwd("/u/awilso10/MOD35")
|
|
| 18 |
|
|
| 14 | 19 |
### directory containing daily files |
| 15 | 20 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 16 | 21 |
|
| ... | ... | |
| 38 | 43 |
fdly$year=format(fdly$date,"%Y") |
| 39 | 44 |
nrow(fdly) |
| 40 | 45 |
|
| 41 |
## check validity (via npar and ntime) of nc files |
|
| 42 |
#for(i in 1:nrow(fdly)){
|
|
| 43 |
# fdly$ntime[i]<-as.numeric(system(paste("cdo -s ntime ",fdly$path[i]),intern=T))
|
|
| 44 |
# fdly$npar[i]<-as.numeric(system(paste("cdo -s npar ",fdly$path[i]),intern=T))
|
|
| 45 |
# fdly$fyear[i]<-as.numeric(system(paste("cdo -s showyear ",fdly$path[i]),intern=T))
|
|
| 46 |
# fdly$fmonth[i]<-as.numeric(system(paste("cdo -s showmon ",fdly$path[i]),intern=T))
|
|
| 47 |
# fdly$fvar[i]<-system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
|
|
| 48 |
# print(paste(i," out of ",nrow(fdly)," for year ", fdly$fyear[i])) |
|
| 49 |
#} |
|
| 50 |
|
|
| 51 | 46 |
## print some summaries |
| 52 | 47 |
if(verbose) print("Summary of available daily files")
|
| 53 | 48 |
print(table(fdly$year)) |
| 54 | 49 |
print(table(fdly$month)) |
| 55 | 50 |
#print(table(fdly$fvar)) |
| 56 | 51 |
|
| 57 |
## Identify which files failed test |
|
| 58 |
#fdly$drop=is.na(fdly$npar)|fdly$fvar!=finalvars |
|
| 59 |
|
|
| 60 |
## delete files that fail check? |
|
| 61 |
delete=F |
|
| 62 |
if(delete) {
|
|
| 63 |
print(paste(sum(fdly$drop),"files will be deleted")) |
|
| 64 |
file.remove(as.character(fdly$path[fdly$drop])) |
|
| 65 |
} |
|
| 66 |
## remove dropped files from list |
|
| 67 |
#fdly=fdly[!fdly$drop,] |
|
| 68 |
|
|
| 69 | 52 |
################################################################################# |
| 70 | 53 |
## Combine the year-by-year files into a single daily file in the summary directory (for archiving) |
| 71 | 54 |
if(verbose) print("Merging daily files into single file output")
|
| ... | ... | |
| 108 | 91 |
|
| 109 | 92 |
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15) |
| 110 | 93 |
|
| 111 |
## subset dates |
|
| 112 |
## due to bug (?) in CDO tools, only 10 years of data can be processed at a time or strange areas of NAs appear. |
|
| 113 |
datesubset="-seldate,2000-01-01,2011-12-31" |
|
| 114 |
|
|
| 115 | 94 |
## Monthly means |
| 116 | 95 |
if(verbose) print("Calculating the monthly means")
|
| 117 | 96 |
system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
| ... | ... | |
| 127 | 106 |
#months=c("01","02","03","04","05","06","07","08","09","10","11","12")
|
| 128 | 107 |
# month="02" |
| 129 | 108 |
|
| 130 |
#ymonmean=function(month){
|
|
| 131 |
# tfile=paste(tempdir(),"/",tile,"_",month,"_",Sys.getpid(),".txt",sep="") |
|
| 132 |
# write.table(paste(fdly$path[fdly$month==month],collapse=" "),tfile,col.names=F,row.names=F,quote=F) |
|
| 133 |
# system(paste("cat ",tfile," | ",ncopath,"ncra -O -o ",tsdir,"/MOD35_",tile,"_",month,".nc",sep=""))
|
|
| 134 |
# system(paste("cdo -O -setyear,",myear," -setmon,",month," -setday,15 -mulc,-1 -subc,100 ",tsdir,"/MOD35_",tile,"_",month,".nc ",tsdir,"/MOD35_",tile,"_",month,"b.nc",sep=""))
|
|
| 135 |
# system(paste(ncopath,"ncrename -v PClear,PCloud ",tsdir,"/MOD35_",tile,"_",month,"b.nc",sep="")) |
|
| 136 |
# system(paste(ncopath,"ncatted ", |
|
| 137 |
#" -a long_name,PCloud,o,c,\"Mean Probability of Cloud\" ", |
|
| 138 |
#" -a missing_value,PCloud,o,b,255 ", |
|
| 139 |
#" -a _FillValue,PCloud,d,b,255 ", |
|
| 140 |
#tsdir,"/MOD35_",tile,"_",month,"b.nc",sep="")) |
|
| 141 |
#} |
|
| 142 |
#mclapply(months,ymonmean) |
|
| 143 |
|
|
| 144 |
## merge to a single file |
|
| 145 |
# system(paste("cdo -O -b I8 -v -mergetime ",paste(tsdir,"/MOD35_",tile,"_",months,"b.nc ",sep="",collapse=" ")," ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
| 146 |
|
|
| 147 |
#system(paste("cdo -v -O selindexbox,1,100,1100,1200 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_dailysmall.nc",sep=""),wait=T)
|
|
| 148 |
#system(paste("cdo -v -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_dailysmall.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
| 149 |
#system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -timmean -selmon,2 ",datesubset," ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
| 150 |
#system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""),wait=T)
|
|
| 151 |
#system(paste("cdo -O sorttimestamp -setyear,",myear," -mulc,-1 -subc,100 -ydrunmean,30 ",datesubset," ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ydrunmean30.nc",sep=""),wait=T)
|
|
| 152 |
#system(paste("scp ",tsdir,"/MOD35_",tile,"_ymonmean.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))
|
|
| 109 |
#system(paste("cdo -O sorttimestamp -setyear,",myear," -mulc,-1 -subc,100 -ydrunmean,30 ",outdir2,"/MOD35_",tile,"_daily.nc ",outdir2,"/MOD35_",tile,"_ydrunmean30.nc &",sep=""),wait=T)
|
|
| 153 | 110 |
#system(paste("scp summary/MOD35_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))
|
| 154 | 111 |
#system(paste("ncdump -h ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",sep=""))
|
| 155 | 112 |
|
| 156 | 113 |
|
| 157 | 114 |
## Monthly standard deviation |
| 158 | 115 |
if(verbose) print("Calculating the monthly SD")
|
| 159 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
| 116 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -monmean ",
|
|
| 117 |
outdir2,"/MOD35_",tile,"_daily.nc ", |
|
| 118 |
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep="")) |
|
| 160 | 119 |
system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep="")) |
| 161 | 120 |
system(paste(ncopath,"ncatted ", |
| 162 | 121 |
" -a long_name,PCloud_sd,o,c,\"Standard Deviation of p(cloud)\" ", |
Also available in: Unified diff
updated MOD35 processing script to work in separate temporary directories and transfer daily cloud data to lou for archiving in preparation for global processing (which cannot be done in personal directory). Also set up MOD35_Climatology script to be submitted to LDAN queue to run on lou rather than Pleiades. However environment still not set up correctly and script fails when reading netcdf files. Andrew is working on it.