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