|
1
|
### 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
|
|
