|
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
|
## report on daily file:
|
|
90
|
#ncfile=paste(outdir2,"/MOD35_",tile,"_daily.nc",sep="")
|
|
91
|
system(paste("ncdump -h ",outdir2,"/MOD35_",tile,"_daily.nc | head -20 ",sep=""))
|
|
92
|
system(paste("cdo showyear ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
93
|
system(paste("cdo showmon ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
94
|
#system(paste("cdo showdate ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
95
|
|
|
96
|
|
|
97
|
### produce a monthly timeseries?
|
|
98
|
#system(paste("cdo -O monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""))
|
|
99
|
|
|
100
|
#############################
|
|
101
|
## Generate the Climatologies
|
|
102
|
if(verbose) print("Generate monthly climatologies")
|
|
103
|
|
|
104
|
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15)
|
|
105
|
|
|
106
|
## subset dates
|
|
107
|
## due to bug (?) in CDO tools, only 10 years of data can be processed at a time or strange areas of NAs appear.
|
|
108
|
datesubset="-seldate,2002-01-01,2012-01-01"
|
|
109
|
|
|
110
|
## Monthly means
|
|
111
|
if(verbose) print("Calculating the monthly means")
|
|
112
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
113
|
|
|
114
|
#system(paste("cdo -v -O selindexbox,1,100,1100,1200 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_dailysmall.nc",sep=""),wait=T)
|
|
115
|
#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)
|
|
116
|
|
|
117
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
118
|
|
|
119
|
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)
|
|
120
|
|
|
121
|
|
|
122
|
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)
|
|
123
|
|
|
124
|
system(paste(ncopath,"ncrename -v PClear,PCloud ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
125
|
system(paste(ncopath,"ncatted ",
|
|
126
|
" -a long_name,PCloud,o,c,\"Mean Probability of Cloud\" ",
|
|
127
|
tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
128
|
|
|
129
|
## Monthly standard deviation
|
|
130
|
if(verbose) print("Calculating the monthly SD")
|
|
131
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
132
|
system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
133
|
system(paste(ncopath,"ncatted ",
|
|
134
|
" -a long_name,PCloud_sd,o,c,\"Standard Deviation of p(cloud)\" ",
|
|
135
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
136
|
|
|
137
|
## frequency of cloud days p(clear<90%)
|
|
138
|
if(verbose) print("Calculating the proportion of cloudy and probably cloudy days")
|
|
139
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -lec,90 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
140
|
system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
141
|
system(paste(ncopath,"ncatted ",
|
|
142
|
" -a long_name,CF,o,c,\"Cloud Frequency: Proportion of Days with probability of clear < 90%\" ",
|
|
143
|
" -a units,CF,o,c,\"Proportion (%)\" ",
|
|
144
|
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
145
|
|
|
146
|
## number of observations
|
|
147
|
if(verbose) print("Calculating the number of missing variables")
|
|
148
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -eqc,9999 -setmisstoc,9999 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
149
|
system(paste(ncopath,"ncrename -v PClear,Pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
150
|
system(paste(ncopath,"ncatted ",
|
|
151
|
" -a long_name,Pmiss,o,c,\"Proportion of Days with missing data\" ",
|
|
152
|
" -a units,Pmiss,o,c,\"Proportion (%)\" ",
|
|
153
|
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
154
|
|
|
155
|
## TODO: fix projection information so GDAL can read it correctly.
|
|
156
|
## clean up variables?
|
|
157
|
|
|
158
|
## append variables to a single file
|
|
159
|
if(verbose) print("Append all monthly climatologies into a single file")
|
|
160
|
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
161
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
162
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
163
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
164
|
|
|
165
|
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
|
|
166
|
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")
|
|
167
|
|
|
168
|
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep=""))
|
|
169
|
#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=""))
|
|
170
|
|
|
171
|
## invert latitude so it plays nicely with gdal
|
|
172
|
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))
|
|
173
|
|
|
174
|
## proj string taken from GDAL-written MODIS tile
|
|
175
|
projstring="PROJCS[\"Sinusoidal (Sanson-Flamsteed)\",GEOGCS[\"wgs84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Sinusoidal\"],PARAMETER[\"longitude_of_center\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]"
|
|
176
|
|
|
177
|
## update attributes
|
|
178
|
system(paste(ncopath,"ncatted ",
|
|
179
|
" -a false_easting,sinusoidal,o,d,0. ",
|
|
180
|
" -a false_northing,sinusoidal,o,d,0. ",
|
|
181
|
" -a longitude_of_central_meridian,sinusoidal,o,d,0. ",
|
|
182
|
" -a longitude_of_prime_meridian,sinusoidal,o,d,0. ",
|
|
183
|
" -a semi_major_axis,sinusoidal,o,d,6378137. ",
|
|
184
|
" -a inverse_flattening,sinusoidal,o,d,298.257223563 ",
|
|
185
|
" -a spatial_ref,sinusoidal,o,c,",projstring,
|
|
186
|
" -a GeoTransform,sinusoidal,o,c,\"-7783653.638366 926.6254331391661 0 1111950.519767 0 -926.6254331391667\" ",
|
|
187
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
|
188
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",
|
|
189
|
" -a tile,global,o,c,\"",tile,"\" ",
|
|
190
|
" -a institution,global,o,c,\"Yale University\" ",
|
|
191
|
" -a source,global,o,c,\"MODIS Cloud Product (MOD35) Collection 6\" ",
|
|
192
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
|
193
|
outdir2,"/MOD35_",tile,".nc",sep=""))
|
|
194
|
|
|
195
|
|
|
196
|
print(paste("############################### Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))
|
|
197
|
print("Years:")
|
|
198
|
print(table(fdly$fyear))
|
|
199
|
print("Months:")
|
|
200
|
print(table(fdly$fmonth))
|
|
201
|
|
|
202
|
## quit R
|
|
203
|
q("no")
|
|
204
|
|
