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
|
require(RSQLite)
|
21
|
|
22
|
## get options
|
23
|
opta <- getopt(matrix(c(
|
24
|
'date', 'd', 1, 'character',
|
25
|
'tile', 't', 1, 'character',
|
26
|
'verbose','v',1,'logical',
|
27
|
'help', 'h', 0, 'logical'
|
28
|
), ncol=4, byrow=TRUE))
|
29
|
if ( !is.null(opta$help) )
|
30
|
{
|
31
|
prg <- commandArgs()[1];
|
32
|
cat(paste("Usage: ", prg, " --date | -d <file> :: The date to process\n", sep=""));
|
33
|
q(status=1);
|
34
|
}
|
35
|
|
36
|
|
37
|
## default date and tile to play with (will be overwritten below when running in batch)
|
38
|
date="20000410"
|
39
|
tile="h11v08"
|
40
|
platform="pleiades"
|
41
|
verbose=T
|
42
|
|
43
|
## now update using options if given
|
44
|
date=opta$date
|
45
|
tile=opta$tile
|
46
|
verbose=opta$verbose #print out extensive information for debugging?
|
47
|
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
48
|
## get year and doy from date
|
49
|
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
50
|
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
51
|
|
52
|
if(platform=="pleiades"){
|
53
|
## location of MOD06 files
|
54
|
datadir=paste("/nobackupp4/datapool/modis/MOD35_L2.006/",year,"/",doy,"/",sep="")
|
55
|
## path to some executables
|
56
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
57
|
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif"
|
58
|
## path to swath database
|
59
|
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"
|
60
|
## specify working directory
|
61
|
setwd("/nobackupp1/awilso10/mod35")
|
62
|
gisBase="/u/armichae/pr/grass-6.4.2/"
|
63
|
## path to MOD11A1 file for this tile to align grid/extent
|
64
|
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
65
|
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
66
|
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 "
|
67
|
}
|
68
|
|
69
|
if(platform=="litoria"){ #if running on local server, use different paths
|
70
|
## specify working directory
|
71
|
setwd("~/acrobates/adamw/projects/interp")
|
72
|
gisBase="/usr/lib/grass64"
|
73
|
## location of MOD06 files
|
74
|
datadir="~/acrobates/adamw/projects/interp/data/modis/mod35"
|
75
|
## path to some executables
|
76
|
ncopath=""
|
77
|
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
|
78
|
## path to swath database
|
79
|
db="~/acrobates/adamw/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
|
80
|
## get grid file
|
81
|
td=raster(paste("~/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
82
|
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 "
|
83
|
}
|
84
|
|
85
|
|
86
|
### print some status messages
|
87
|
if(verbose) print(paste("Processing tile",tile," for date",date))
|
88
|
|
89
|
## load tile information and get bounding box
|
90
|
load(file="modlandTiles.Rdata")
|
91
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
92
|
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner
|
93
|
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner
|
94
|
|
95
|
|
96
|
## vars to process
|
97
|
vars=as.data.frame(matrix(c(
|
98
|
"Cloud_Mask", "CM",
|
99
|
"Quality_Assurance", "QA"),
|
100
|
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
101
|
|
102
|
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end.
|
103
|
finalvars=c("PClear")
|
104
|
|
105
|
|
106
|
#####################################################
|
107
|
##find swaths in region from sqlite database for the specified date/tile
|
108
|
if(verbose) print("Accessing swath ID's from database")
|
109
|
con=dbConnect("SQLite", dbname = db)
|
110
|
fs=dbGetQuery(con,paste("SELECT * from swath_geo6
|
111
|
WHERE east>=",tile_bb$lon_min," AND
|
112
|
west<=",tile_bb$lon_max," AND
|
113
|
north>=",tile_bb$lat_min," AND
|
114
|
south<=",tile_bb$lat_max," AND
|
115
|
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND
|
116
|
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j"))
|
117
|
))
|
118
|
con=dbDisconnect(con)
|
119
|
fs$id=substr(fs$id,7,19)
|
120
|
## find the swaths on disk (using datadir)
|
121
|
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T)
|
122
|
|
123
|
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk"))
|
124
|
|
125
|
|
126
|
############################################################################
|
127
|
############################################################################
|
128
|
### Use the HEG tool to grid all available swath data for this date-tile
|
129
|
for(file in swaths){
|
130
|
## Function to generate hegtool parameter file for multi-band HDF-EOS file
|
131
|
print(paste("Starting file",basename(file)))
|
132
|
outfile=paste(tempdir(),"/",basename(file),sep="")
|
133
|
## First write the parameter file (careful, heg is very finicky!)
|
134
|
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
135
|
grp=paste("
|
136
|
BEGIN
|
137
|
INPUT_FILENAME=",file,"
|
138
|
OBJECT_NAME=mod35
|
139
|
FIELD_NAME=",vars$variable,"|
|
140
|
BAND_NUMBER = 1
|
141
|
OUTPUT_PIXEL_SIZE_X=1000
|
142
|
OUTPUT_PIXEL_SIZE_Y=1000
|
143
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
144
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
145
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
146
|
RESAMPLING_TYPE =NN
|
147
|
OUTPUT_PROJECTION_TYPE = SIN
|
148
|
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
|
149
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
150
|
ELLIPSOID_CODE = WGS84
|
151
|
OUTPUT_TYPE = HDFEOS
|
152
|
OUTPUT_FILENAME = ",outfile,"
|
153
|
END
|
154
|
|
155
|
",sep="")
|
156
|
|
157
|
## if any remnants from previous runs remain, delete them
|
158
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
159
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
160
|
## write it to a file
|
161
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
162
|
## now run the swath2grid tool
|
163
|
## write the gridded file
|
164
|
log=system(paste("(",swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
165
|
## clean up temporary files in working directory
|
166
|
# file.remove(list.files(pattern=
|
167
|
# paste("filetable.temp_",
|
168
|
# as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process
|
169
|
if(verbose) print(log)
|
170
|
print(paste("Finished gridding ", file))
|
171
|
} #end looping over swaths
|
172
|
|
173
|
########################
|
174
|
## confirm at least one file for this date is present. If not, quit.
|
175
|
outfiles=paste(tempdir(),"/",basename(swaths),sep="")
|
176
|
if(!any(file.exists(outfiles))) {
|
177
|
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
178
|
q("no",status=0)
|
179
|
}
|
180
|
|
181
|
#####################################################
|
182
|
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths
|
183
|
|
184
|
## Identify output file
|
185
|
ncfile=paste(outdir,"/MOD35_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
186
|
|
187
|
## function to convert binary to decimal to assist in identifying correct values
|
188
|
## this is helpful when defining QA handling below, but isn't used in processing
|
189
|
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
190
|
## for example:
|
191
|
## b2d(c(T,T))
|
192
|
|
193
|
## set Grass to overwrite
|
194
|
Sys.setenv(GRASS_OVERWRITE=1)
|
195
|
Sys.setenv(DEBUG=1)
|
196
|
Sys.setenv(GRASS_GUI="txt")
|
197
|
|
198
|
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
199
|
## make temporary working directory
|
200
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
201
|
if(!file.exists(tf)) dir.create(tf)
|
202
|
## create output directory if needed
|
203
|
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
|
204
|
|
205
|
## set up temporary grass instance for this PID
|
206
|
if(verbose) print(paste("Set up temporary grass session in",tf))
|
207
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
208
|
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
209
|
|
210
|
## Define region by importing one MOD11A1 raster.
|
211
|
print("Import one MOD11A1 raster to define grid")
|
212
|
if(platform=="pleiades") {
|
213
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
214
|
output="modisgrid",flags=c("quiet","overwrite","o"))
|
215
|
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
216
|
}
|
217
|
|
218
|
if(platform=="litoria"){
|
219
|
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
220
|
output="modisgrid",flags=c("overwrite","o"))
|
221
|
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
222
|
}
|
223
|
## Identify which files to process
|
224
|
tfs=basename(swaths)
|
225
|
## drop swaths that did not produce an output file (typically due to not overlapping the ROI)
|
226
|
tfs=tfs[tfs%in%list.files(tempdir())]
|
227
|
#tfs=list.files(tempdir(),pattern="temp.*hdf")
|
228
|
nfs=length(tfs)
|
229
|
if(verbose) print(paste(nfs,"swaths available for processing"))
|
230
|
|
231
|
## loop through scenes and process QA flags
|
232
|
for(i in 1:nfs){
|
233
|
file=paste(tempdir(),"/",tfs[i],sep="")
|
234
|
## Cloud Mask
|
235
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
|
236
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
237
|
## QA ## extract first bit to keep only "useful" values of cloud mask
|
238
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
|
239
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
240
|
## extract first two bits of cloud mask
|
241
|
system(paste("r.mapcalc <<EOF
|
242
|
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0)
|
243
|
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999)
|
244
|
Pclear1_",i," = if(CM_cloud_",i,"==0,0,if(CM_cloud_",i,"==1,66,if(CM_cloud_",i,"==2,95,if(CM_cloud_",i,"==3,99,-9999))))
|
245
|
Pclear_",i," = if(QA_useful_",i,"==1,Pclear1_",i,",-9999)
|
246
|
EOF",sep=""))
|
247
|
|
248
|
# CM_path_",i," = ((CM1_",i," / 2^6) % 2^2)
|
249
|
|
250
|
execGRASS("r.null",map=paste("Pclear_",i,sep=""),setnull="-9999")
|
251
|
execGRASS("r.null",map=paste("CM_cloud_",i,sep=""),setnull="-9999")
|
252
|
|
253
|
|
254
|
} #end loop through sub daily files
|
255
|
|
256
|
#### Now generate daily minimum p(clear)
|
257
|
system(paste("r.mapcalc <<EOF
|
258
|
Pclear_daily=int((min(",paste("if(isnull(Pclear_",1:nfs,"),9999,Pclear_",1:nfs,")",sep="",collapse=","),")))
|
259
|
EOF",sep=""))
|
260
|
|
261
|
|
262
|
## reset null values
|
263
|
execGRASS("r.null",map="Pclear_daily",setnull="9999")
|
264
|
|
265
|
|
266
|
### Write the files to a netcdf file
|
267
|
## create image group to facilitate export as multiband netcdf
|
268
|
execGRASS("i.group",group="mod35",input=c("Pclear_daily")) ; print("")
|
269
|
|
270
|
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
271
|
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
272
|
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
273
|
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
274
|
|
275
|
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
276
|
## create temporary nc file with time information to append to MOD06 data
|
277
|
cat(paste("
|
278
|
netcdf time {
|
279
|
dimensions:
|
280
|
time = 1 ;
|
281
|
variables:
|
282
|
int time(time) ;
|
283
|
time:units = \"days since 2000-01-01 00:00:00\" ;
|
284
|
time:calendar = \"gregorian\";
|
285
|
time:long_name = \"time of observation\";
|
286
|
data:
|
287
|
time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";
|
288
|
}"),file=paste(tempdir(),"/time.cdl",sep=""))
|
289
|
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
290
|
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
|
291
|
## add other attributes
|
292
|
system(paste(ncopath,"ncrename -v Band1,PClear ",ncfile,sep=""))
|
293
|
system(paste(ncopath,"ncatted -a scale_factor,PClear,o,d,1 -a units,PClear,o,c,\"Probability (%)\" -a missing_value,PClear,o,d,255 -a _FillValue,PClear,o,d,255 -a long_name,PClear,o,c,\"Probability of Clear Sky\" ",ncfile,sep=""))
|
294
|
|
295
|
# system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
|
296
|
|
297
|
|
298
|
### delete the temporary files
|
299
|
unlink_.gislock()
|
300
|
system(paste("rm -frR ",tf,sep=""))
|
301
|
|
302
|
|
303
|
## Confirm that the file has the correct attributes, otherwise delete it
|
304
|
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
305
|
## confirm it has all 'final variables as specified above"
|
306
|
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
|
307
|
|
308
|
if(ntime!=1|!fvar) {
|
309
|
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
310
|
file.remove(ncfile)
|
311
|
}
|
312
|
|
313
|
## print out some info
|
314
|
print(paste(" ################################################################### Finished ",date,
|
315
|
"################################################################"))
|
316
|
|
317
|
## delete old files
|
318
|
#system("cleartemp")
|
319
|
|
320
|
## quit
|
321
|
q("no",status=0)
|