1
|
###################################################################################
|
2
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
3
|
|
4
|
|
5
|
## connect to server of choice
|
6
|
#system("ssh litoria")
|
7
|
#R
|
8
|
|
9
|
library(sp)
|
10
|
library(rgdal)
|
11
|
library(reshape)
|
12
|
library(ncdf4)
|
13
|
library(geosphere)
|
14
|
library(rgeos)
|
15
|
library(multicore)
|
16
|
library(raster)
|
17
|
library(lattice)
|
18
|
library(rgl)
|
19
|
library(hdf5)
|
20
|
library(spgrass6)
|
21
|
|
22
|
X11.options(type="Xlib")
|
23
|
ncores=20 #number of threads to use
|
24
|
|
25
|
setwd("/home/adamw/personal/projects/interp")
|
26
|
setwd("/home/adamw/acrobates/projects/interp")
|
27
|
|
28
|
roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")
|
29
|
roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
30
|
|
31
|
### Downloading data from LAADSWEB
|
32
|
# subset by geographic area of interest
|
33
|
# subset: 40-47, -115--125
|
34
|
|
35
|
## download data from ftp site. Unfortunately the selection has to be selected using the website and orders downloaded via ftp.
|
36
|
|
37
|
system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")
|
38
|
|
39
|
## specify some working directories
|
40
|
gdir="output/"
|
41
|
datadir="data/modis/MOD06_L2_hdf"
|
42
|
outdir="data/modis/MOD06_L2_hdf2"
|
43
|
tifdir="/media/data/MOD06_L2_tif"
|
44
|
summarydatadir="data/modis/MOD06_climatologies"
|
45
|
|
46
|
|
47
|
fs=data.frame(
|
48
|
path=list.files(datadir,full=T,recursive=T,pattern="hdf"),
|
49
|
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))
|
50
|
fs$date=as.Date(substr(fs$file,11,17),"%Y%j")
|
51
|
fs$month=format(fs$date,"%m")
|
52
|
fs$year=format(fs$date,"%Y")
|
53
|
fs$time=substr(fs$file,19,22)
|
54
|
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))
|
55
|
fs$dateid=format(fs$date,"%Y%m%d")
|
56
|
fs$path=as.character(fs$path)
|
57
|
fs$file=as.character(fs$file)
|
58
|
|
59
|
## output ROI
|
60
|
#get bounding box of region in m
|
61
|
ge=SpatialPoints(data.frame(lon=c(-125,-115),lat=c(40,47)))
|
62
|
projection(ge)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
63
|
ge2=spTransform(ge, CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
64
|
|
65
|
|
66
|
## vars
|
67
|
vars=as.data.frame(matrix(c(
|
68
|
"Cloud_Effective_Radius", "CER",
|
69
|
"Cloud_Effective_Radius_Uncertainty", "CERU",
|
70
|
"Cloud_Optical_Thickness", "COT",
|
71
|
"Cloud_Optical_Thickness_Uncertainty", "COTU",
|
72
|
"Cloud_Water_Path", "CWP",
|
73
|
"Cloud_Water_Path_Uncertainty", "CWPU",
|
74
|
"Cloud_Phase_Optical_Properties", "CPOP",
|
75
|
"Cloud_Multi_Layer_Flag", "CMLF",
|
76
|
"Cloud_Mask_1km", "CM1",
|
77
|
"Quality_Assurance_1km", "QA"),
|
78
|
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)
|
79
|
|
80
|
|
81
|
### Installation of hegtool
|
82
|
## needed 32-bit libraries and java for program to install correctly
|
83
|
|
84
|
# system(paste("hegtool -h ",fs$path[1],sep=""))
|
85
|
|
86
|
|
87
|
#### Function to generate hegtool parameter file for multi-band HDF-EOS file
|
88
|
swath2grid=function(i=1,files,vars,outdir,upleft="47 -125",lowright="41 -115"){
|
89
|
file=fs$path[i]
|
90
|
print(paste("Starting file",basename(file)))
|
91
|
outfile=paste(outdir,"/",fs$file[i],sep="")
|
92
|
# date=fs$date[1]
|
93
|
# origin=as.POSIXct("1970-01-01 00:00:00",tz="GMT")
|
94
|
### First write the parameter file (careful, heg is very finicky!)
|
95
|
hdr=paste("NUM_RUNS = ",length(vars),"|MULTI_BAND_HDFEOS:",length(vars),sep="")
|
96
|
grp=paste("
|
97
|
BEGIN
|
98
|
INPUT_FILENAME=",file,"
|
99
|
OBJECT_NAME=mod06
|
100
|
FIELD_NAME=",vars,"|
|
101
|
BAND_NUMBER = 1
|
102
|
OUTPUT_PIXEL_SIZE_X=1000
|
103
|
OUTPUT_PIXEL_SIZE_Y=1000
|
104
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
105
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
106
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
107
|
RESAMPLING_TYPE =NN
|
108
|
OUTPUT_PROJECTION_TYPE = SIN
|
109
|
ELLIPSOID_CODE = WGS84
|
110
|
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 86400.0 0.0 0.0 )
|
111
|
OUTPUT_TYPE = HDFEOS
|
112
|
OUTPUT_FILENAME = ",outfile,"
|
113
|
END
|
114
|
|
115
|
",sep="")
|
116
|
## if any remnants from previous runs remain, delete them
|
117
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
118
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
119
|
## write it to a file
|
120
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
121
|
## now run the swath2grid tool - must be run as root (argh!)!
|
122
|
## write the tiff file
|
123
|
log=system(paste("sudo MRTDATADIR=/usr/local/heg/data ",
|
124
|
"PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",
|
125
|
paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
126
|
print(paste("Finished ", file))
|
127
|
}
|
128
|
|
129
|
|
130
|
### update fs with completed files
|
131
|
fs$complete=fs$file%in%list.files(outdir,pattern="hdf$")
|
132
|
table(fs$complete)
|
133
|
|
134
|
#### Run the gridding procedure
|
135
|
|
136
|
system("sudo ls")
|
137
|
|
138
|
mclapply(which(!fs$complete),function(fi){
|
139
|
swath2grid(fi,vars=vars$variable,files=fs,
|
140
|
outdir=outdir,
|
141
|
upleft="47 -125",lowright="40 -115")},
|
142
|
mc.cores=24)
|
143
|
|
144
|
|
145
|
##############################################################
|
146
|
### Import to GRASS for processing
|
147
|
|
148
|
#fs$grass=paste(fs$month,fs$year,fs$file,sep="_")
|
149
|
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""))
|
150
|
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 "
|
151
|
|
152
|
## fucntion to convert binary to decimal to assist in identifying correct values
|
153
|
b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
154
|
## for example:
|
155
|
b2d(c(T,T))
|
156
|
|
157
|
### create (or connect to) grass location
|
158
|
gisDbase="/media/data/grassdata"
|
159
|
gisLocation="oregon"
|
160
|
gisMapset="mod06"
|
161
|
## set Grass to overwrite
|
162
|
Sys.setenv(GRASS_OVERWRITE=1)
|
163
|
Sys.setenv(DEBUG=0)
|
164
|
|
165
|
## temporary objects to test function below
|
166
|
i=1
|
167
|
file=paste(outdir,"/",fs$file[1],sep="")
|
168
|
date=as.Date("2000-03-02")
|
169
|
|
170
|
|
171
|
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
172
|
loadcloud<-function(date,fs){
|
173
|
### set up grass session
|
174
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")
|
175
|
|
176
|
## set up tempfile for this PID
|
177
|
initGRASS(gisBase="/usr/lib/grass64",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
178
|
system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +datum=WGS84 +units=m +no_defs\"",sep=""))
|
179
|
|
180
|
#system("g.mapset PERMANENT")
|
181
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),
|
182
|
output="modisgrid",flags=c("quiet","overwrite","o"))
|
183
|
system("g.region rast=modisgrid save=roi --overwrite")
|
184
|
system("g.region roi")
|
185
|
system("g.region -p")
|
186
|
# getLocationProj()
|
187
|
|
188
|
|
189
|
## Identify which files to process
|
190
|
tfs=fs$file[fs$date==date]
|
191
|
nfs=length(tfs)
|
192
|
|
193
|
### print some summary info
|
194
|
print(date)
|
195
|
## loop through scenes and process QA flags
|
196
|
for(i in 1:nfs){
|
197
|
file=paste(outdir,"/",tfs[i],sep="")
|
198
|
## Cloud Mask
|
199
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
200
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
201
|
## extract cloudy and 'confidently clear' pixels
|
202
|
system(paste("r.mapcalc <<EOF
|
203
|
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0
|
204
|
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3
|
205
|
EOF",sep=""))
|
206
|
|
207
|
## QA
|
208
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
209
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
210
|
## QA_CER
|
211
|
system(paste("r.mapcalc <<EOF
|
212
|
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
213
|
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3
|
214
|
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0
|
215
|
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1
|
216
|
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3
|
217
|
EOF",sep=""))
|
218
|
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1
|
219
|
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3
|
220
|
|
221
|
## Optical Thickness
|
222
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
223
|
output=paste("COT_",i,sep=""),
|
224
|
title="cloud_effective_radius",
|
225
|
flags=c("overwrite","o")) ; print("")
|
226
|
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
227
|
## keep only positive COT values where quality is 'useful' and 'very good' & scale to real units
|
228
|
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,"*0.009999999776482582,null())\"",sep=""))
|
229
|
## set COT to 0 in clear-sky pixels
|
230
|
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
231
|
|
232
|
## Effective radius ##
|
233
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
234
|
output=paste("CER_",i,sep=""),
|
235
|
title="cloud_effective_radius",
|
236
|
flags=c("overwrite","o")) ; print("")
|
237
|
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
238
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
239
|
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
240
|
## set CER to 0 in clear-sky pixels
|
241
|
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
242
|
|
243
|
## Cloud Water Path
|
244
|
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
245
|
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
246
|
# flags=c("overwrite","o")) ; print("")
|
247
|
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
248
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
249
|
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))
|
250
|
## set CER to 0 in clear-sky pixels
|
251
|
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
252
|
|
253
|
|
254
|
} #end loop through sub daily files
|
255
|
|
256
|
#### Now generate daily averages (or maximum in case of cloud flag)
|
257
|
|
258
|
system(paste("r.mapcalc <<EOF
|
259
|
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
260
|
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
261
|
COT_daily=COT_numer/COT_denom
|
262
|
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
263
|
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
264
|
CER_daily=CER_numer/CER_denom
|
265
|
CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),")
|
266
|
EOF",sep=""))
|
267
|
|
268
|
#### Write the file to a geotiff
|
269
|
execGRASS("r.out.gdal",input="CER_daily",output=paste(tifdir,"/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))
|
270
|
execGRASS("r.out.gdal",input="COT_daily",output=paste(tifdir,"/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))
|
271
|
execGRASS("r.out.gdal",input="CLD_daily",output=paste(tifdir,"/CLD_",format(date,"%Y%m%d"),".tif",sep=""),nodata=99,flags=c("quiet"))
|
272
|
|
273
|
### delete the temporary files
|
274
|
unlink_.gislock()
|
275
|
system("/usr/lib/grass64/etc/clean_temp")
|
276
|
system(paste("rm -R ",tf,sep=""))
|
277
|
### print update
|
278
|
print(paste(" ################################################################### Finished ",date,"
|
279
|
################################################################"))
|
280
|
}
|
281
|
|
282
|
|
283
|
###########################################
|
284
|
### Now run it
|
285
|
|
286
|
tdates=sort(unique(fs$date))
|
287
|
done=tdates%in%as.Date(substr(list.files(tifdir),5,12),"%Y%m%d")
|
288
|
table(done)
|
289
|
tdates=tdates[!done]
|
290
|
|
291
|
mclapply(tdates,function(date) loadcloud(date,fs=fs))
|
292
|
|
293
|
|
294
|
|
295
|
#######################################################################################33
|
296
|
### Produce the monthly averages
|
297
|
|
298
|
## get list of daily files
|
299
|
fs2=data.frame(
|
300
|
path=list.files(tifdir,full=T,recursive=T,pattern="tif$"),
|
301
|
file=basename(list.files(tifdir,full=F,recursive=T,pattern="tif$")))
|
302
|
fs2$type=substr(fs2$file,1,3)
|
303
|
fs2$date=as.Date(substr(fs2$file,5,12),"%Y%m%d")
|
304
|
fs2$month=format(fs2$date,"%m")
|
305
|
fs2$year=format(fs2$date,"%Y")
|
306
|
fs2$path=as.character(fs2$path)
|
307
|
fs2$file=as.character(fs2$file)
|
308
|
|
309
|
|
310
|
# Define type/month products
|
311
|
vs=expand.grid(type=unique(fs2$type),month=c("01","02","03","04","05","06","07","08","09","10","11","12"))
|
312
|
|
313
|
## identify which have been completed
|
314
|
#done=
|
315
|
# do.call(rbind,strsplit(list.files(summarydatadir),"_|[.]"))[,3]
|
316
|
#table(done)
|
317
|
#tdates=tdates[!done]
|
318
|
|
319
|
|
320
|
## process the summaries using the raster package
|
321
|
mclapply(1:nrow(vs),function(i){
|
322
|
print(paste("Starting ",vs$type[i]," for month ",vs$month[i]))
|
323
|
td=stack(fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])])
|
324
|
print(paste("Processing Metric ",vs$type[i]," for month ",vs$month[i]))
|
325
|
calc(td,mean,na.rm=T,
|
326
|
filename=paste(summarydatadir,"/",vs$type[i],"_mean_",vs$month[i],".tif",sep=""),
|
327
|
format="GTiff")
|
328
|
calc(td,sd,na.rm=T,
|
329
|
filename=paste(summarydatadir,"/",vs$type[i],"_sd_",vs$month[i],".tif",sep=""),
|
330
|
format="GTiff")
|
331
|
print(paste("Processing missing data for ",vs$type[i]," for month ",vs$month[i]))
|
332
|
calc(td,function(i)
|
333
|
sum(!is.na(i)),filename=paste(summarydatadir,"/",vs$type[i],"_count_",vs$month[i],".tif",sep=""),
|
334
|
format="GTiff")
|
335
|
calc(td,function(i) sum(ifelse(i==0,0,1)),
|
336
|
filename=paste(summarydatadir,"/",vs$type[i],"_clear_",vs$month[i],".tif",sep=""),format="GTiff")
|
337
|
gc()
|
338
|
}
|
339
|
)
|
340
|
|
341
|
|