|
1
|
### Process a folder of daily MOD35 HDF files to produce a climatology
|
|
2
|
|
|
3
|
.libPaths("/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
|
4
|
|
|
5
|
## import commandline arguments
|
|
6
|
library(getopt,lib="/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
|
7
|
## get options
|
|
8
|
opta <- getopt(matrix(c(
|
|
9
|
'tile', 't', 1, 'character',
|
|
10
|
'verbose','v',1,'logical'
|
|
11
|
), ncol=4, byrow=TRUE))
|
|
12
|
|
|
13
|
tile=opta$tile #tile="h00v08"
|
|
14
|
verbose=opta$verbose #print out extensive information for debugging?
|
|
15
|
|
|
16
|
## set working directory
|
|
17
|
setwd("/nobackupp1/awilso10/mod35")
|
|
18
|
#setwd("/u/awilso10/MOD35")
|
|
19
|
|
|
20
|
### directory containing daily files
|
|
21
|
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
|
22
|
|
|
23
|
### directory to hold climatology
|
|
24
|
outdir2="summary" #directory for combined daily files and summarized files
|
|
25
|
if(!file.exists(outdir2)) dir.create(outdir2)
|
|
26
|
|
|
27
|
### path to NCO
|
|
28
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
|
29
|
|
|
30
|
################################################################################
|
|
31
|
## Get list of all daily files
|
|
32
|
if(verbose) print(paste("Checking daily output in preparation for generating climatology:",tile))
|
|
33
|
|
|
34
|
fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F)
|
|
35
|
fdly$file=basename(fdly$path)
|
|
36
|
fdly$dateid=substr(fdly$file,14,21)
|
|
37
|
fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")
|
|
38
|
fdly$month=format(fdly$date,"%m")
|
|
39
|
fdly$year=format(fdly$date,"%Y")
|
|
40
|
nrow(fdly)
|
|
41
|
|
|
42
|
## print some summaries
|
|
43
|
if(verbose) print(paste("Summary of available daily files:",tile))
|
|
44
|
print(table(fdly$year))
|
|
45
|
print(table(fdly$month))
|
|
46
|
#print(table(fdly$fvar))
|
|
47
|
|
|
48
|
#################################################################################
|
|
49
|
## Combine the year-by-year files into a single daily file in the summary directory (for archiving)
|
|
50
|
if(verbose) print(paste("Merging daily files into single file output:",tile))
|
|
51
|
|
|
52
|
## create temporary directory to put intermediate files (will be deleted when R quits)
|
|
53
|
tsdir=paste(tempdir(),"/summary",sep="")
|
|
54
|
if(!file.exists(tsdir)) dir.create(tsdir,recursive=T)
|
|
55
|
|
|
56
|
## merge all daily files to create a single file with all dates
|
|
57
|
system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
58
|
system(paste("ncdump -h ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
59
|
|
|
60
|
## Update attributes
|
|
61
|
system(paste(ncopath,"ncatted ",
|
|
62
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Summaries\" ",
|
|
63
|
" -a institution,global,o,c,\"Yale University\" ",
|
|
64
|
" -a source,global,o,c,\"MODIS Collection 6 Cloud Mask (MOD35)\" ",
|
|
65
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
|
66
|
outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
67
|
|
|
68
|
## report on daily file:
|
|
69
|
#ncfile=paste(outdir2,"/MOD35_",tile,"_daily.nc",sep="")
|
|
70
|
system(paste("ncdump -h ",outdir2,"/MOD35_",tile,"_daily.nc | head -20 ",sep=""))
|
|
71
|
system(paste("cdo showyear ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
72
|
system(paste("cdo showmon ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
73
|
#system(paste("cdo showdate ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
|
74
|
|
|
75
|
|
|
76
|
### produce a monthly timeseries?
|
|
77
|
#system(paste("cdo -O monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""))
|
|
78
|
|
|
79
|
#############################
|
|
80
|
## Generate the Climatologies
|
|
81
|
if(verbose) print(paste("Generate monthly climatologies: ",tile))
|
|
82
|
|
|
83
|
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15)
|
|
84
|
|
|
85
|
## Overall Means
|
|
86
|
if(verbose) print(paste("Calculating the overall mean:",tile))
|
|
87
|
system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setmon,1 -setday,1 -mulc,100 -timmean -lec,1 ",outdir2,"/MOD35_",tile,"_daily.nc ",outdir2,"/MOD35_",tile,"_mean.nc",sep=""),wait=T)
|
|
88
|
system(paste(ncopath,"ncrename -v CMday,CFday -v CMnight,CFnight ",outdir2,"/MOD35_",tile,"_mean.nc",sep=""))
|
|
89
|
system(paste(ncopath,"ncatted ",
|
|
90
|
" -a long_name,CFday,o,c,\"Daytime Cloud Frequency\" ",
|
|
91
|
" -a missing_value,CFday,o,b,255 ",
|
|
92
|
" -a _FillValue,CFday,d,b,255 ",
|
|
93
|
" -a long_name,CFnight,o,c,\"Nighttime Cloud Frequency\" ",
|
|
94
|
" -a missing_value,CFnight,o,b,255 ",
|
|
95
|
" -a _FillValue,CFnight,d,b,255 ",
|
|
96
|
outdir2,"/MOD35_",tile,"_mean.nc",sep=""))
|
|
97
|
|
|
98
|
## Monthly means
|
|
99
|
if(verbose) print(paste("Calculating the monthly means:",tile))
|
|
100
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -mulc,100 -ymonmean -lec,1 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
101
|
system(paste(ncopath,"ncrename -v CMday,CFday -v CMnight,CFnight ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
102
|
system(paste(ncopath,"ncatted ",
|
|
103
|
" -a long_name,CFday,o,c,\"Daytime Cloud Frequency\" ",
|
|
104
|
" -a units,CFday,o,c,\"Proportion (%)\" ",
|
|
105
|
" -a missing_value,CFday,o,b,255 ",
|
|
106
|
" -a _FillValue,CFday,d,b,255 ",
|
|
107
|
" -a long_name,CFnight,o,c,\"Nighttime Cloud Frequency\" ",
|
|
108
|
" -a units,CFnight,o,c,\"Proportion (%)\" ",
|
|
109
|
" -a missing_value,CFnight,o,b,255 ",
|
|
110
|
" -a _FillValue,CFnight,d,b,255 ",
|
|
111
|
tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
112
|
|
|
113
|
|
|
114
|
## Monthly Mean
|
|
115
|
#months=c("01","02","03","04","05","06","07","08","09","10","11","12")
|
|
116
|
# month="02"
|
|
117
|
|
|
118
|
#system(paste("cdo -O sorttimestamp -setyear,",myear," -mulc,-1 -subc,100 -ydrunmean,30 ",outdir2,"/MOD35_",tile,"_daily.nc ",outdir2,"/MOD35_",tile,"_ydrunmean30.nc &",sep=""),wait=T)
|
|
119
|
#system(paste("scp summary/MOD35_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))
|
|
120
|
#system(paste("ncdump -h ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",sep=""))
|
|
121
|
|
|
122
|
## Monthly standard deviation
|
|
123
|
if(verbose) print(paste("Calculating the monthly SD:",tile))
|
|
124
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -mulc,100 -monmean -lec,1 ",
|
|
125
|
outdir2,"/MOD35_",tile,"_daily.nc ",
|
|
126
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
127
|
system(paste(ncopath,"ncrename -v CMday,CFday_sd -v CMnight,CFnight_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
128
|
system(paste(ncopath,"ncatted ",
|
|
129
|
" -a long_name,CFnight_sd,o,c,\"Standard Deviation of monthly nighttime cloud frequency\" ",
|
|
130
|
" -a long_name,CFday_sd,o,c,\"Standard Deviation of monthly daytime cloud frequency\" ",
|
|
131
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
132
|
|
|
133
|
## frequency of cloud days p(clear<90%)
|
|
134
|
#if(verbose) print(paste("Calculating the proportion of cloudy and probably cloudy days:",tile))
|
|
135
|
#system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -lec,90 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
136
|
#system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
137
|
#system(paste(ncopath,"ncatted ",
|
|
138
|
#" -a long_name,CF,o,c,\"Cloud Frequency: Proportion of Days with probability of clear < 90%\" ",
|
|
139
|
#" -a units,CF,o,c,\"Proportion (%)\" ",
|
|
140
|
#tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
141
|
|
|
142
|
## number of observations
|
|
143
|
if(verbose) print(paste("Calculating the number of missing variables:",tile))
|
|
144
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -eqc,9999 -setmisstoc,9999 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
145
|
system(paste(ncopath,"ncrename -v CMday,CFday_pmiss -v CMnight,CFnight_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
146
|
system(paste(ncopath,"ncatted ",
|
|
147
|
" -a long_name,CFday_pmiss,o,c,\"Proportion of Days with missing data\" ",
|
|
148
|
" -a units,CFday_pmiss,o,c,\"Proportion (%)\" ",
|
|
149
|
" -a long_name,CFnight_pmiss,o,c,\"Proportion of Days with missing data\" ",
|
|
150
|
" -a units,CFnight_pmiss,o,c,\"Proportion (%)\" ",
|
|
151
|
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
152
|
|
|
153
|
## TODO: fix projection information so GDAL can read it correctly.
|
|
154
|
## clean up variables?
|
|
155
|
|
|
156
|
## append variables to a single file
|
|
157
|
if(verbose) print(paste("Append all monthly climatologies into a single file:",tile))
|
|
158
|
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
159
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
160
|
#system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
161
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
162
|
|
|
163
|
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
|
|
164
|
if(verbose) print(paste("Clean up file (update attributes, flip latitudes, add grid description:",tile))
|
|
165
|
|
|
166
|
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep=""))
|
|
167
|
#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=""))
|
|
168
|
|
|
169
|
## invert latitude so it plays nicely with gdal
|
|
170
|
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))
|
|
171
|
|
|
172
|
## proj string taken from GDAL-written MODIS tile
|
|
173
|
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]]"
|
|
174
|
|
|
175
|
## update attributes
|
|
176
|
system(paste(ncopath,"ncatted ",
|
|
177
|
" -a false_easting,sinusoidal,o,d,0. ",
|
|
178
|
" -a false_northing,sinusoidal,o,d,0. ",
|
|
179
|
" -a longitude_of_central_meridian,sinusoidal,o,d,0. ",
|
|
180
|
" -a longitude_of_prime_meridian,sinusoidal,o,d,0. ",
|
|
181
|
" -a semi_major_axis,sinusoidal,o,d,6378137. ",
|
|
182
|
" -a inverse_flattening,sinusoidal,o,d,298.257223563 ",
|
|
183
|
" -a spatial_ref,sinusoidal,o,c,",projstring,
|
|
184
|
" -a GeoTransform,sinusoidal,o,c,\"-7783653.638366 926.6254331391661 0 1111950.519767 0 -926.6254331391667\" ",
|
|
185
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
|
186
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",
|
|
187
|
" -a tile,global,o,c,\"",tile,"\" ",
|
|
188
|
" -a institution,global,o,c,\"Yale University\" ",
|
|
189
|
" -a source,global,o,c,\"MODIS Cloud Product (MOD35) Collection 6\" ",
|
|
190
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
|
191
|
outdir2,"/MOD35_",tile,".nc",sep=""))
|
|
192
|
|
|
193
|
|
|
194
|
print(paste("############################### Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))
|
|
195
|
print("Years:")
|
|
196
|
print(table(fdly$fyear))
|
|
197
|
print("Months:")
|
|
198
|
print(table(fdly$fmonth))
|
|
199
|
|
|
200
|
## quit R
|
|
201
|
q("no")
|
|
202
|
|
