Revision f27b378d
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD06_Climatology.r | ||
|---|---|---|
| 13 | 13 |
|
| 14 | 14 |
### directory containing daily files |
| 15 | 15 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 16 |
if(!file.exists(outdir)) dir.create(outdir) |
|
| 17 | 16 |
|
| 18 | 17 |
### directory to hold climatology |
| 19 | 18 |
outdir2="summary" #directory for combined daily files and summarized files |
| ... | ... | |
| 31 | 30 |
## Get list of all daily files |
| 32 | 31 |
if(verbose) print("Checking daily output in preparation for generating climatology")
|
| 33 | 32 |
|
| 34 |
fdly=data.frame( |
|
| 35 |
path=list.files(outdir,pattern="nc$",full=T), |
|
| 36 |
file=list.files(outdir,pattern="nc$")) |
|
| 33 |
fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F) |
|
| 34 |
fdly$file=basename(fdly$path) |
|
| 37 | 35 |
fdly$dateid=substr(fdly$file,14,21) |
| 38 | 36 |
fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d") |
| 39 | 37 |
fdly$month=format(fdly$date,"%m") |
| ... | ... | |
| 41 | 39 |
nrow(fdly) |
| 42 | 40 |
|
| 43 | 41 |
## check validity (via npar and ntime) of nc files |
| 44 |
for(i in 1:nrow(fdly)){
|
|
| 45 |
fdly$ntime[i]<-as.numeric(system(paste("cdo -s ntime ",fdly$path[i]),intern=T))
|
|
| 46 |
fdly$npar[i]<-as.numeric(system(paste("cdo -s npar ",fdly$path[i]),intern=T))
|
|
| 47 |
fdly$fyear[i]<-as.numeric(system(paste("cdo -s showyear ",fdly$path[i]),intern=T))
|
|
| 48 |
fdly$fmonth[i]<-as.numeric(system(paste("cdo -s showmon ",fdly$path[i]),intern=T))
|
|
| 49 |
fdly$fvar[i]<-system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
|
|
| 50 |
print(paste(i," out of ",nrow(fdly)," for year ", fdly$fyear[i])) |
|
| 51 |
} |
|
| 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 |
#}
|
|
| 52 | 50 |
|
| 53 | 51 |
## print some summaries |
| 54 | 52 |
if(verbose) print("Summary of available daily files")
|
| 55 |
print(table(fdly$fyear))
|
|
| 56 |
print(table(fdly$fmonth))
|
|
| 57 |
print(table(fdly$fvar)) |
|
| 53 |
print(table(fdly$year)) |
|
| 54 |
print(table(fdly$month)) |
|
| 55 |
#print(table(fdly$fvar))
|
|
| 58 | 56 |
|
| 59 | 57 |
## Identify which files failed test |
| 60 |
fdly$drop=is.na(fdly$npar)|fdly$fvar!=finalvars |
|
| 58 |
#fdly$drop=is.na(fdly$npar)|fdly$fvar!=finalvars
|
|
| 61 | 59 |
|
| 62 | 60 |
## delete files that fail check? |
| 63 |
delete=T
|
|
| 61 |
delete=F
|
|
| 64 | 62 |
if(delete) {
|
| 65 | 63 |
print(paste(sum(fdly$drop),"files will be deleted")) |
| 66 | 64 |
file.remove(as.character(fdly$path[fdly$drop])) |
| 67 | 65 |
} |
| 68 | 66 |
## remove dropped files from list |
| 69 |
fdly=fdly[!fdly$drop,] |
|
| 67 |
#fdly=fdly[!fdly$drop,]
|
|
| 70 | 68 |
|
| 71 | 69 |
################################################################################# |
| 72 | 70 |
## Combine the year-by-year files into a single daily file in the summary directory (for archiving) |
| ... | ... | |
| 98 | 96 |
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15) |
| 99 | 97 |
|
| 100 | 98 |
## Monthly means |
| 99 |
if(verbose) print("Calculating the monthly means")
|
|
| 101 | 100 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""),wait=T)
|
| 102 | 101 |
|
| 103 | 102 |
## Monthly standard deviation |
| 103 |
if(verbose) print("Calculating the monthly SD")
|
|
| 104 | 104 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
|
| 105 | 105 |
system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep="")) |
| 106 | 106 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 110 | 110 |
tsdir,"/MOD06_",tile,"_ymonstd.nc",sep="")) |
| 111 | 111 |
|
| 112 | 112 |
## cer > 20 |
| 113 |
if(verbose) print("Calculating the proportion of days with CER > 20")
|
|
| 113 | 114 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -gtc,20 -selvar,CER ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
|
| 114 | 115 |
system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep="")) |
| 115 | 116 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 117 | 118 |
" -a units,CER20,o,c,\"Proportion\" ", |
| 118 | 119 |
tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep="")) |
| 119 | 120 |
|
| 121 |
## cld == 0 |
|
| 122 |
if(verbose) print("Calculating the proportion of cloudy days")
|
|
| 123 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,0 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld0.nc",sep=""))
|
|
| 124 |
system(paste(ncopath,"ncrename -v CLD2,CLD0 ",tsdir,"/MOD06_",tile,"_ymoncld0.nc",sep="")) |
|
| 125 |
system(paste(ncopath,"ncatted ", |
|
| 126 |
" -a long_name,CLD0,o,c,\"Proportion of Days with Cloud Mask == 0\" ", |
|
| 127 |
" -a units,CLD0,o,c,\"Proportion\" ", |
|
| 128 |
tsdir,"/MOD06_",tile,"_ymoncld0.nc",sep="")) |
|
| 129 |
|
|
| 130 |
## cld == 0|1 |
|
| 131 |
if(verbose) print("Calculating the proportion of cloudy days")
|
|
| 132 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -lec,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld01.nc",sep=""))
|
|
| 133 |
system(paste(ncopath,"ncrename -v CLD2,CLD01 ",tsdir,"/MOD06_",tile,"_ymoncld01.nc",sep="")) |
|
| 134 |
system(paste(ncopath,"ncatted ", |
|
| 135 |
" -a long_name,CLD01,o,c,\"Proportion of Days with Cloud Mask == 0|1\" ", |
|
| 136 |
" -a units,CLD01,o,c,\"Proportion\" ", |
|
| 137 |
tsdir,"/MOD06_",tile,"_ymoncld01.nc",sep="")) |
|
| 138 |
|
|
| 139 |
## cld == 1 |
|
| 140 |
if(verbose) print("Calculating the proportion of uncertain days")
|
|
| 141 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld1.nc",sep=""))
|
|
| 142 |
system(paste(ncopath,"ncrename -v CLD2,CLD1 ",tsdir,"/MOD06_",tile,"_ymoncld1.nc",sep="")) |
|
| 143 |
system(paste(ncopath,"ncatted ", |
|
| 144 |
" -a long_name,CLD1,o,c,\"Proportion of Days with Cloud Mask == 1 (uncertain)\" ", |
|
| 145 |
" -a units,CLD1,o,c,\"Proportion\" ", |
|
| 146 |
tsdir,"/MOD06_",tile,"_ymoncld1.nc",sep="")) |
|
| 147 |
|
|
| 148 |
|
|
| 149 |
## cld >= 2 (setting cld==01 to missing because 'uncertain') |
|
| 150 |
if(verbose) print("Calculating the proportion of clear days")
|
|
| 151 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -gtc,1 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld2.nc",sep=""))
|
|
| 152 |
system(paste(ncopath,"ncrename -v CLD2,CLD23 ",tsdir,"/MOD06_",tile,"_ymoncld2.nc",sep="")) |
|
| 153 |
system(paste(ncopath,"ncatted ", |
|
| 154 |
" -a long_name,CLD23,o,c,\"Proportion of Days with Cloud Mask >= 2 (Probably Clear or Certainly Clear)\" ", |
|
| 155 |
" -a units,CLD23,o,c,\"Proportion\" ", |
|
| 156 |
tsdir,"/MOD06_",tile,"_ymoncld2.nc",sep="")) |
|
| 157 |
|
|
| 158 |
## cld >= 1 |
|
| 159 |
if(verbose) print("Calculating the proportion of clear days")
|
|
| 160 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -gec,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld13.nc",sep=""))
|
|
| 161 |
system(paste(ncopath,"ncrename -v CLD2,CLD13 ",tsdir,"/MOD06_",tile,"_ymoncld13.nc",sep="")) |
|
| 162 |
system(paste(ncopath,"ncatted ", |
|
| 163 |
" -a long_name,CLD13,o,c,\"Proportion of Days with Cloud Mask >= 1\" ", |
|
| 164 |
" -a units,CLD13,o,c,\"Proportion\" ", |
|
| 165 |
tsdir,"/MOD06_",tile,"_ymoncld13.nc",sep="")) |
|
| 166 |
|
|
| 120 | 167 |
## number of observations |
| 168 |
if(verbose) print("Calculating the number of missing variables")
|
|
| 121 | 169 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,9999 -setmisstoc,9999 -selvar,CER,CLD ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
|
| 122 | 170 |
system(paste(ncopath,"ncrename -v CER,CER_pmiss -v CLD,CLD_pmiss ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep="")) |
| 123 | 171 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 134 | 182 |
system(paste(ncopath,"ncks -O ",tsdir,"/MOD06_",tile,"_ymonmean.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
| 135 | 183 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonstd.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
| 136 | 184 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncer20.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
| 185 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncld0.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 186 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncld1.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 187 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncld2.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 137 | 188 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
| 138 | 189 |
|
| 139 | 190 |
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes |
Also available in: Unified diff
Updates to MOD06 Climatology Script and MDO06 processing to isolate the effects of the MOD35 -landcover bias