Project

General

Profile

Download (11.1 KB) Statistics
| Branch: | Revision:
1 f9d40a14 Adam M. Wilson @ pfe
### Process a folder of daily MOD06 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 f27b378d Adam M. Wilson @ pfe
 fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F)
34
  fdly$file=basename(fdly$path)
35 f9d40a14 Adam M. Wilson @ pfe
  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 f27b378d Adam M. Wilson @ pfe
#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 f9d40a14 Adam M. Wilson @ pfe
51
## print some summaries
52
if(verbose) print("Summary of available daily files")
53 f27b378d Adam M. Wilson @ pfe
print(table(fdly$year))
54
print(table(fdly$month))
55
#print(table(fdly$fvar))
56 f9d40a14 Adam M. Wilson @ pfe
57 9c8b4de6 Adam M. Wilson @ pfe
## Identify which files failed test
58 f27b378d Adam M. Wilson @ pfe
#fdly$drop=is.na(fdly$npar)|fdly$fvar!=finalvars
59 9c8b4de6 Adam M. Wilson @ pfe
60 f9d40a14 Adam M. Wilson @ pfe
## delete files that fail check?
61 f27b378d Adam M. Wilson @ pfe
delete=F
62 f9d40a14 Adam M. Wilson @ pfe
if(delete) {
63
  print(paste(sum(fdly$drop),"files will be deleted"))
64
  file.remove(as.character(fdly$path[fdly$drop]))
65
}
66 9c8b4de6 Adam M. Wilson @ pfe
## remove dropped files from list
67 f27b378d Adam M. Wilson @ pfe
#fdly=fdly[!fdly$drop,]
68 f9d40a14 Adam M. Wilson @ pfe
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,"/MOD06_",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 (MOD06) Daily Timeseries\" ",
84
" -a institution,global,o,c,\"Yale University\" ",
85
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ",
86
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
87
outdir2,"/MOD06_",tile,"_daily.nc",sep=""))
88
89
### produce a monthly timeseries?
90
#system(paste("cdo -O monmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",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 f27b378d Adam M. Wilson @ pfe
if(verbose) print("Calculating the monthly means")
100 f9d40a14 Adam M. Wilson @ pfe
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""),wait=T)
101
102
## Monthly standard deviation
103 f27b378d Adam M. Wilson @ pfe
if(verbose) print("Calculating the monthly SD")
104 f9d40a14 Adam M. Wilson @ pfe
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
105
system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",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,"/MOD06_",tile,"_ymonstd.nc",sep=""))
111
112
## cer > 20
113 f27b378d Adam M. Wilson @ pfe
if(verbose) print("Calculating the proportion of days with CER > 20")
114 f9d40a14 Adam M. Wilson @ pfe
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=""))
115
system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
116
system(paste(ncopath,"ncatted ",
117
" -a long_name,CER20,o,c,\"Proportion of Days with Cloud Particle Effective Radius > 20um\" ",
118
" -a units,CER20,o,c,\"Proportion\" ",
119
tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
120
121 f27b378d Adam M. Wilson @ pfe
## 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
167 f9d40a14 Adam M. Wilson @ pfe
## number of observations
168 f27b378d Adam M. Wilson @ pfe
if(verbose) print("Calculating the number of missing variables")
169 f9d40a14 Adam M. Wilson @ pfe
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=""))
170
system(paste(ncopath,"ncrename -v CER,CER_pmiss -v CLD,CLD_pmiss ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
171
system(paste(ncopath,"ncatted ",
172
" -a long_name,CER_pmiss,o,c,\"Proportion of Days with missing data for CER\" ",
173
" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",
174
" -a scale_factor,CER_pmiss,o,d,0.01 ",
175
" -a units,CER_pmiss,o,c,\"Proportion\" ",
176
" -a scale_factor,CLD_pmiss,o,d,0.01 ",
177
" -a units,CLD_pmiss,o,c,\"Proportion\" ",
178
tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
179
180
## append variables to a single file
181
if(verbose) print("Append all monthly climatologies into a single file")
182
system(paste(ncopath,"ncks -O ",tsdir,"/MOD06_",tile,"_ymonmean.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
183
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonstd.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
184
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncer20.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
185 f27b378d Adam M. Wilson @ pfe
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=""))
188 f9d40a14 Adam M. Wilson @ pfe
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
189
190
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
191
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")
192
193
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1],"  ",tsdir,"/sinusoidal.nc",sep=""))
194
#system(paste(ncopath,"ncks -A -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1],"  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
195
196
## invert latitude so it plays nicely with gdal
197
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD06_",tile,"_ymon.nc ",outdir2,"/MOD06_",tile,".nc",sep=""))
198
199
## update attributes
200
system(paste(ncopath,"ncatted ",
201
#" -a standard_parallel,sinusoidal,o,c,\"0\" ",
202
#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
203
#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
204
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
205
" -a title,global,o,c,\"MODIS Cloud Product (MOD06) Climatology\" ",
206
" -a institution,global,o,c,\"Yale University\" ",
207
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ",
208
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
209
outdir2,"/MOD06_",tile,".nc",sep=""))
210
211
212
print(paste("###############################  Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))
213
print("Years:")
214
print(table(fdly$fyear))
215
print("Months:")
216
print(table(fdly$fmonth))
217
218
## quit R
219
q("no")
220