| 1 |
35d59dc1
|
Adam M. Wilson @ pfe
|
### Process a folder of daily MOD35 HDF files to produce a climatology
|
| 2 |
|
|
|
| 3 |
|
|
## import commandline arguments
|
| 4 |
|
|
library(getopt)
|
| 5 |
|
|
## get options
|
| 6 |
|
|
opta <- getopt(matrix(c(
|
| 7 |
|
|
'tile', 't', 1, 'character',
|
| 8 |
|
|
'verbose','v',1,'logical'
|
| 9 |
|
|
), ncol=4, byrow=TRUE))
|
| 10 |
|
|
|
| 11 |
|
|
tile=opta$tile #tile="h11v08"
|
| 12 |
|
|
verbose=opta$verbose #print out extensive information for debugging?
|
| 13 |
|
|
|
| 14 |
|
|
### directory containing daily files
|
| 15 |
|
|
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 16 |
|
|
|
| 17 |
|
|
### directory to hold climatology
|
| 18 |
|
|
outdir2="summary" #directory for combined daily files and summarized files
|
| 19 |
|
|
if(!file.exists(outdir2)) dir.create(outdir2)
|
| 20 |
|
|
|
| 21 |
|
|
### path to NCO
|
| 22 |
|
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
| 23 |
|
|
|
| 24 |
|
|
### Vector of variables that must be in file or they will be deleted.
|
| 25 |
|
|
### Formated as output from system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
|
| 26 |
|
|
#finalvars=" CER COT CLD"
|
| 27 |
|
|
|
| 28 |
|
|
|
| 29 |
|
|
################################################################################
|
| 30 |
|
|
## Get list of all daily files
|
| 31 |
|
|
if(verbose) print("Checking daily output in preparation for generating climatology")
|
| 32 |
|
|
|
| 33 |
|
|
fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F)
|
| 34 |
|
|
fdly$file=basename(fdly$path)
|
| 35 |
|
|
fdly$dateid=substr(fdly$file,14,21)
|
| 36 |
|
|
fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")
|
| 37 |
|
|
fdly$month=format(fdly$date,"%m")
|
| 38 |
|
|
fdly$year=format(fdly$date,"%Y")
|
| 39 |
|
|
nrow(fdly)
|
| 40 |
|
|
|
| 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 |
|
|
## print some summaries
|
| 52 |
|
|
if(verbose) print("Summary of available daily files")
|
| 53 |
|
|
print(table(fdly$year))
|
| 54 |
|
|
print(table(fdly$month))
|
| 55 |
|
|
#print(table(fdly$fvar))
|
| 56 |
|
|
|
| 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 |
|
|
#################################################################################
|
| 70 |
|
|
## Combine the year-by-year files into a single daily file in the summary directory (for archiving)
|
| 71 |
|
|
if(verbose) print("Merging daily files into single file output")
|
| 72 |
|
|
|
| 73 |
|
|
## create temporary directory to put intermediate files (will be deleted when R quits)
|
| 74 |
|
|
tsdir=paste(tempdir(),"/summary",sep="")
|
| 75 |
|
|
if(!file.exists(tsdir)) dir.create(tsdir,recursive=T)
|
| 76 |
|
|
|
| 77 |
|
|
## merge all daily files to create a single file with all dates
|
| 78 |
|
|
system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 79 |
|
|
|
| 80 |
|
|
## Update attributes
|
| 81 |
|
|
system(paste(ncopath,"ncatted ",
|
| 82 |
|
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
| 83 |
|
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Daily Timeseries\" ",
|
| 84 |
|
|
" -a institution,global,o,c,\"Yale University\" ",
|
| 85 |
|
|
" -a source,global,o,c,\"MODIS Cloud Mask (MOD35)\" ",
|
| 86 |
|
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
| 87 |
|
|
outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 88 |
|
|
|
| 89 |
|
|
### produce a monthly timeseries?
|
| 90 |
|
|
#system(paste("cdo -O monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""))
|
| 91 |
|
|
|
| 92 |
|
|
#############################
|
| 93 |
|
|
## Generate the Climatologies
|
| 94 |
|
|
if(verbose) print("Generate monthly climatologies")
|
| 95 |
|
|
|
| 96 |
|
|
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15)
|
| 97 |
|
|
|
| 98 |
|
|
## Monthly means
|
| 99 |
|
|
if(verbose) print("Calculating the monthly means")
|
| 100 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,1.0 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
| 101 |
|
|
|
| 102 |
|
|
## Monthly standard deviation
|
| 103 |
|
|
if(verbose) print("Calculating the monthly SD")
|
| 104 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 105 |
|
|
system(paste(ncopath,"ncrename -v CLD,CLD_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 106 |
|
|
system(paste(ncopath,"ncatted ",
|
| 107 |
|
|
#" -a long_name,CER_sd,o,c,\"Cloud Particle Effective Radius (standard deviation of daily observations)\" ",
|
| 108 |
|
|
" -a long_name,CLD_sd,o,c,\"Cloud Mask (standard deviation of daily observations)\" ",
|
| 109 |
|
|
#" -a long_name,COT_sd,o,c,\"Cloud Optical Thickness (standard deviation of daily observations)\" ",
|
| 110 |
|
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 111 |
|
|
|
| 112 |
|
|
## cld == 0
|
| 113 |
|
|
if(verbose) print("Calculating the proportion of cloudy days")
|
| 114 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,0 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))
|
| 115 |
|
|
system(paste(ncopath,"ncrename -v CLD2,CLD0 ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))
|
| 116 |
|
|
system(paste(ncopath,"ncatted ",
|
| 117 |
|
|
" -a long_name,CLD0,o,c,\"Proportion of Days with Cloud Mask == 0\" ",
|
| 118 |
|
|
" -a units,CLD0,o,c,\"Proportion\" ",
|
| 119 |
|
|
tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))
|
| 120 |
|
|
|
| 121 |
|
|
## cld == 0|1
|
| 122 |
|
|
if(verbose) print("Calculating the proportion of cloudy days")
|
| 123 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -lec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 124 |
|
|
system(paste(ncopath,"ncrename -v CLD2,CLD01 ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 125 |
|
|
system(paste(ncopath,"ncatted ",
|
| 126 |
|
|
" -a long_name,CLD01,o,c,\"Proportion of Days with Cloud Mask == 0|1\" ",
|
| 127 |
|
|
" -a units,CLD01,o,c,\"Proportion\" ",
|
| 128 |
|
|
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 129 |
|
|
|
| 130 |
|
|
## cld == 1
|
| 131 |
|
|
if(verbose) print("Calculating the proportion of uncertain days")
|
| 132 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))
|
| 133 |
|
|
system(paste(ncopath,"ncrename -v CLD2,CLD1 ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))
|
| 134 |
|
|
system(paste(ncopath,"ncatted ",
|
| 135 |
|
|
" -a long_name,CLD1,o,c,\"Proportion of Days with Cloud Mask == 1 (uncertain)\" ",
|
| 136 |
|
|
" -a units,CLD1,o,c,\"Proportion\" ",
|
| 137 |
|
|
tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))
|
| 138 |
|
|
|
| 139 |
|
|
|
| 140 |
|
|
## cld >= 2 (setting cld==01 to missing because 'uncertain')
|
| 141 |
|
|
if(verbose) print("Calculating the proportion of clear days")
|
| 142 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -gtc,1 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))
|
| 143 |
|
|
system(paste(ncopath,"ncrename -v CLD2,CLD23 ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))
|
| 144 |
|
|
system(paste(ncopath,"ncatted ",
|
| 145 |
|
|
" -a long_name,CLD23,o,c,\"Proportion of Days with Cloud Mask >= 2 (Probably Clear or Certainly Clear)\" ",
|
| 146 |
|
|
" -a units,CLD23,o,c,\"Proportion\" ",
|
| 147 |
|
|
tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))
|
| 148 |
|
|
|
| 149 |
|
|
## cld >= 1
|
| 150 |
|
|
if(verbose) print("Calculating the proportion of clear days")
|
| 151 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -gec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))
|
| 152 |
|
|
system(paste(ncopath,"ncrename -v CLD2,CLD13 ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))
|
| 153 |
|
|
system(paste(ncopath,"ncatted ",
|
| 154 |
|
|
" -a long_name,CLD13,o,c,\"Proportion of Days with Cloud Mask >= 1\" ",
|
| 155 |
|
|
" -a units,CLD13,o,c,\"Proportion\" ",
|
| 156 |
|
|
tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))
|
| 157 |
|
|
|
| 158 |
|
|
## number of observations
|
| 159 |
|
|
if(verbose) print("Calculating the number of missing variables")
|
| 160 |
|
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,9999 -setmisstoc,9999 -selvar,CLD ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 161 |
|
|
system(paste(ncopath,"ncrename -v CLD,CLD_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 162 |
|
|
system(paste(ncopath,"ncatted ",
|
| 163 |
|
|
" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",
|
| 164 |
|
|
" -a scale_factor,CLD_pmiss,o,d,0.01 ",
|
| 165 |
|
|
" -a units,CLD_pmiss,o,c,\"Proportion\" ",
|
| 166 |
|
|
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 167 |
|
|
|
| 168 |
|
|
## TODO: fix projection information so GDAL can read it correctly.
|
| 169 |
|
|
## clean up variables?
|
| 170 |
|
|
|
| 171 |
|
|
## append variables to a single file
|
| 172 |
|
|
if(verbose) print("Append all monthly climatologies into a single file")
|
| 173 |
|
|
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 174 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 175 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld0.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 176 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 177 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld1.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 178 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld2.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 179 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 180 |
|
|
|
| 181 |
|
|
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
|
| 182 |
|
|
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")
|
| 183 |
|
|
|
| 184 |
|
|
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep=""))
|
| 185 |
|
|
#system(paste(ncopath,"ncks -A -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 186 |
|
|
|
| 187 |
|
|
## invert latitude so it plays nicely with gdal
|
| 188 |
|
|
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 189 |
|
|
|
| 190 |
|
|
## update attributes
|
| 191 |
|
|
system(paste(ncopath,"ncatted ",
|
| 192 |
|
|
#" -a standard_parallel,sinusoidal,o,c,\"0\" ",
|
| 193 |
|
|
#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
|
| 194 |
|
|
#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
|
| 195 |
|
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
| 196 |
|
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",
|
| 197 |
|
|
" -a institution,global,o,c,\"Yale University\" ",
|
| 198 |
|
|
" -a source,global,o,c,\"MODIS Cloud Product (MOD35)\" ",
|
| 199 |
|
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
| 200 |
|
|
outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 201 |
|
|
|
| 202 |
|
|
|
| 203 |
|
|
print(paste("############################### Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))
|
| 204 |
|
|
print("Years:")
|
| 205 |
|
|
print(table(fdly$fyear))
|
| 206 |
|
|
print("Months:")
|
| 207 |
|
|
print(table(fdly$fmonth))
|
| 208 |
|
|
|
| 209 |
|
|
## quit R
|
| 210 |
|
|
q("no")
|
| 211 |
|
|
|
