|
1
|
###################################################################################
|
|
2
|
### R code to aquire and process MOD35_L2 cloud data from the MODIS platform
|
|
3
|
|
|
4
|
|
|
5
|
# Redirect all warnings to stderr()
|
|
6
|
#options(warn = -1)
|
|
7
|
#write("2) write() to stderr", stderr())
|
|
8
|
#write("2) write() to stdout", stdout())
|
|
9
|
#warning("2) warning()")
|
|
10
|
|
|
11
|
|
|
12
|
## import commandline arguments
|
|
13
|
library(getopt)
|
|
14
|
## load libraries
|
|
15
|
require(reshape)
|
|
16
|
require(geosphere)
|
|
17
|
require(raster)
|
|
18
|
require(rgdal)
|
|
19
|
require(spgrass6)
|
|
20
|
|
|
21
|
## get options
|
|
22
|
opta <- getopt(matrix(c(
|
|
23
|
'date', 'd', 1, 'character',
|
|
24
|
'tile', 't', 1, 'character',
|
|
25
|
'verbose','v',1,'logical',
|
|
26
|
'help', 'h', 0, 'logical'
|
|
27
|
), ncol=4, byrow=TRUE))
|
|
28
|
if ( !is.null(opta$help) )
|
|
29
|
{
|
|
30
|
prg <- commandArgs()[1];
|
|
31
|
cat(paste("Usage: ", prg, " --date | -d <file> :: The date to process\n", sep=""));
|
|
32
|
q(status=1);
|
|
33
|
}
|
|
34
|
|
|
35
|
testing=T
|
|
36
|
platform="pleiades"
|
|
37
|
|
|
38
|
## default date and tile to play with (will be overwritten below when running in batch)
|
|
39
|
if(testing){
|
|
40
|
date="20090129"
|
|
41
|
tile="h11v08"
|
|
42
|
tile="h17v00"
|
|
43
|
verbose=T
|
|
44
|
}
|
|
45
|
|
|
46
|
## now update using options if given
|
|
47
|
if(!testing){
|
|
48
|
date=opta$date
|
|
49
|
tile=opta$tile
|
|
50
|
verbose=opta$verbose #print out extensive information for debugging?
|
|
51
|
}
|
|
52
|
## get year and doy from date
|
|
53
|
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
|
54
|
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
|
55
|
|
|
56
|
if(platform=="pleiades"){
|
|
57
|
## location of MOD35 files
|
|
58
|
datadir=paste("/nobackupp4/datapool/modis/MOD09GA.005/",year,"/",doy,"/",sep="")
|
|
59
|
## path to some executables
|
|
60
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
|
61
|
## specify working directory
|
|
62
|
outdir=paste("/nobackupp1/awilso10/mod09/daily/",tile,"/",sep="") #directory for separate daily files
|
|
63
|
basedir="/nobackupp1/awilso10/mod09/" #directory to hold files temporarily before transferring to lou
|
|
64
|
setwd(tempdir())
|
|
65
|
## grass database
|
|
66
|
gisBase="/u/armichae/pr/grass-6.4.2/"
|
|
67
|
## path to MOD11A1 file for this tile to align grid/extent
|
|
68
|
gridfile=list.files("/nobackupp1/awilso10/mod35/MODTILES/",pattern=tile,full=T)[1]
|
|
69
|
td=raster(gridfile)
|
|
70
|
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "
|
|
71
|
}
|
|
72
|
|
|
73
|
### print some status messages
|
|
74
|
if(verbose) print(paste("Processing tile",tile," for date",date))
|
|
75
|
|
|
76
|
## load tile information and get bounding box
|
|
77
|
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata")
|
|
78
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
|
79
|
|
|
80
|
#####################################################
|
|
81
|
|
|
82
|
## function to convert binary to decimal to assist in identifying correct values
|
|
83
|
## this is helpful when defining QA handling below, but isn't used in processing
|
|
84
|
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
|
85
|
## for example:
|
|
86
|
## b2d(c(T,T))
|
|
87
|
|
|
88
|
## set Grass to overwrite
|
|
89
|
Sys.setenv(GRASS_OVERWRITE=1)
|
|
90
|
Sys.setenv(DEBUG=1)
|
|
91
|
Sys.setenv(GRASS_GUI="txt")
|
|
92
|
|
|
93
|
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
|
94
|
## make temporary working directory
|
|
95
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
|
96
|
if(!file.exists(tf)) dir.create(tf)
|
|
97
|
## create output directory if needed
|
|
98
|
## Identify output file
|
|
99
|
ncfile=paste(outdir,"MOD09cloud_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
|
100
|
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
|
|
101
|
|
|
102
|
## set up temporary grass instance for this PID
|
|
103
|
if(verbose) print(paste("Set up temporary grass session in",tf))
|
|
104
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod09",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
105
|
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
|
106
|
|
|
107
|
## Define region by importing one MOD11A1 raster.
|
|
108
|
print("Import one MOD11A1 raster to define grid")
|
|
109
|
if(platform=="pleiades") {
|
|
110
|
execGRASS("r.in.gdal",input=td@file@name,output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
111
|
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
112
|
}
|
|
113
|
|
|
114
|
if(platform=="litoria"){
|
|
115
|
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
|
116
|
output="modisgrid",flags=c("overwrite","o"))
|
|
117
|
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
118
|
}
|
|
119
|
|
|
120
|
## Identify which files to process
|
|
121
|
tfs=unique(sub("CM_|QA_|SenZen_|SolZen_","",basename(outfiles)))
|
|
122
|
#tfs=list.files(tempdir(),pattern="temp.*hdf")
|
|
123
|
nfs=length(tfs)
|
|
124
|
if(verbose) print(paste(nfs,"swaths available for processing"))
|
|
125
|
|
|
126
|
## loop through scenes and process QA flags
|
|
127
|
for(i in 1:nfs){
|
|
128
|
bfile=tfs[i]
|
|
129
|
## Read in the data from the HDFs
|
|
130
|
## Cloud Mask
|
|
131
|
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
|
132
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
133
|
## QA ## extract first bit to keep only "useful" values of cloud mask
|
|
134
|
execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),
|
|
135
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
136
|
## Sensor Zenith ## extract first bit to keep only "low angle" observations
|
|
137
|
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
|
138
|
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
139
|
|
|
140
|
## check for interpolation artefacts
|
|
141
|
# execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))
|
|
142
|
# execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))
|
|
143
|
# execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))
|
|
144
|
|
|
145
|
## write out the table of weights for use in the neighborhood analysis to identify bad pixels from swtif
|
|
146
|
p=75 #must be odd
|
|
147
|
mat=matrix(rep(0,p*p),nrow=p)
|
|
148
|
mat[0.5+p/2,]=1
|
|
149
|
cat(mat,file="weights.txt")
|
|
150
|
execGRASS("r.neighbors",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"clump",sep=""),method="range",size=p,weight="weights.txt") # too slow!
|
|
151
|
system(paste("r.mapcalc \"SZ_",i,"_clump2=SZ_",i,"clump==0\"",sep=""))
|
|
152
|
|
|
153
|
# p=-50:50
|
|
154
|
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
155
|
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",min(p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
156
|
# vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))
|
|
157
|
# colnames(vals)=c("value","count")
|
|
158
|
# vals$count=as.numeric(as.character(vals$count))
|
|
159
|
# vals$value=as.numeric(as.character(vals$value))
|
|
160
|
# vals=na.omit(vals)
|
|
161
|
# vals$count[vals$value==1&vals$count>10]
|
|
162
|
#
|
|
163
|
#plot(p~value,data=vals)
|
|
164
|
# print(sum(vals$p[vals$p>.1]))
|
|
165
|
|
|
166
|
## Solar Zenith ## extract first bit to keep only "low angle" observations
|
|
167
|
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
|
168
|
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
169
|
## produce the summaries
|
|
170
|
system(paste("r.mapcalc <<EOF
|
|
171
|
CM_fill_",i," = if(isnull(CM1_",i,"),1,0)
|
|
172
|
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0)
|
|
173
|
SZ_low_",i," = if(SZ_",i,"_clump2==0&SZ_",i,"<6000,1,0)
|
|
174
|
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0)
|
|
175
|
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
|
176
|
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null())
|
|
177
|
SZday_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==1,SZ_",i,",null())
|
|
178
|
SZnight_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==0,SZ_",i,",null())
|
|
179
|
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())
|
|
180
|
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())
|
|
181
|
EOF",sep=""))
|
|
182
|
|
|
183
|
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
|
184
|
# CM_dscore_",i," = if((CM_dayflag_",i,"==0|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,if(SoZ_",i,">=8500,3,4))))
|
|
185
|
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))
|
|
186
|
|
|
187
|
drawplot=F
|
|
188
|
if(drawplot){
|
|
189
|
d2=stack(
|
|
190
|
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
|
191
|
raster(readRAST6(paste("CM1_",i,sep=""))),
|
|
192
|
# raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
193
|
# raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
194
|
# raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
195
|
# raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
196
|
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
|
197
|
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
|
198
|
raster(readRAST6(paste("SZ_",i,sep=""))),
|
|
199
|
raster(readRAST6(paste("SZ_",i,"_clump",sep=""))),
|
|
200
|
raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))
|
|
201
|
)
|
|
202
|
plot(d2,add=F)
|
|
203
|
}
|
|
204
|
|
|
205
|
|
|
206
|
} #end loop through sub daily files
|
|
207
|
|
|
208
|
## select lowest view angle
|
|
209
|
## use r.series to find minimum
|
|
210
|
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
|
211
|
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
|
212
|
|
|
213
|
## select cloud observation with lowest sensor zenith for day and night
|
|
214
|
system(
|
|
215
|
paste("r.mapcalc <<EOF
|
|
216
|
CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
|
217
|
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
|
218
|
))
|
|
219
|
|
|
220
|
if(plot){
|
|
221
|
ps=1:nfs
|
|
222
|
ps=c(12,14,17)
|
|
223
|
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))
|
|
224
|
sz_clump=brick(lapply(ps,function(i) raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))))
|
|
225
|
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
226
|
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
|
227
|
library(rasterVis)
|
|
228
|
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))
|
|
229
|
levelplot(sz_clump)
|
|
230
|
levelplot(d)
|
|
231
|
levelplot(d2)
|
|
232
|
}
|
|
233
|
|
|
234
|
|
|
235
|
### Write the files to a netcdf file
|
|
236
|
## create image group to facilitate export as multiband netcdf
|
|
237
|
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily")) ; print("")
|
|
238
|
|
|
239
|
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
|
240
|
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
|
241
|
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
|
242
|
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
|
243
|
|
|
244
|
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
|
245
|
## create temporary nc file with time information to append to MOD06 data
|
|
246
|
cat(paste("
|
|
247
|
netcdf time {
|
|
248
|
dimensions:
|
|
249
|
time = 1 ;
|
|
250
|
variables:
|
|
251
|
int time(time) ;
|
|
252
|
time:units = \"days since 2000-01-01 00:00:00\" ;
|
|
253
|
time:calendar = \"gregorian\";
|
|
254
|
time:long_name = \"time of observation\";
|
|
255
|
data:
|
|
256
|
time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";
|
|
257
|
}"),file=paste(tempdir(),"/time.cdl",sep=""))
|
|
258
|
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
|
259
|
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
|
|
260
|
## add other attributes
|
|
261
|
system(paste(ncopath,"ncrename -v Band1,CMday -v Band2,CMnight ",ncfile,sep=""))
|
|
262
|
system(paste(ncopath,"ncatted ",
|
|
263
|
" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ",
|
|
264
|
" -a missing_value,CMday,o,b,255 ",
|
|
265
|
" -a _FillValue,CMday,o,b,255 ",
|
|
266
|
" -a valid_range,CMday,o,b,\"0,3\" ",
|
|
267
|
" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ",
|
|
268
|
" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ",
|
|
269
|
" -a missing_value,CMnight,o,b,255 ",
|
|
270
|
" -a _FillValue,CMnight,o,b,255 ",
|
|
271
|
" -a valid_range,CMnight,o,b,\"0,3\" ",
|
|
272
|
" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ",
|
|
273
|
ncfile,sep=""))
|
|
274
|
#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
|
|
275
|
|
|
276
|
|
|
277
|
## Confirm that the file has the correct attributes, otherwise delete it
|
|
278
|
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
|
279
|
## confirm it has all 'final variables as specified above"
|
|
280
|
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
|
|
281
|
|
|
282
|
if(ntime!=1|!fvar) {
|
|
283
|
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
|
284
|
file.remove(ncfile)
|
|
285
|
}
|
|
286
|
############ copy files to lou
|
|
287
|
#if(platform=="pleiades"){
|
|
288
|
# archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
|
289
|
# system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
|
290
|
# system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
|
291
|
# file.remove(ncfile)
|
|
292
|
# file.remove(paste(ncfile,".aux.xml",sep=""))
|
|
293
|
#}
|
|
294
|
|
|
295
|
|
|
296
|
### delete the temporary files
|
|
297
|
# unlink_.gislock()
|
|
298
|
# system(paste("rm -frR ",tempdir(),sep=""))
|
|
299
|
|
|
300
|
|
|
301
|
## print out some info
|
|
302
|
print(paste("####################### Finished ",tile,"-",date, "###################################"))
|
|
303
|
|
|
304
|
## delete old files
|
|
305
|
#system("cleartemp")
|
|
306
|
|
|
307
|
## quit
|
|
308
|
q("no",status=0)
|
