|
1
|
###################################################################################
|
|
2
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
|
3
|
|
|
4
|
|
|
5
|
##First read in the arguments listed at the command line
|
|
6
|
args=(commandArgs(TRUE))
|
|
7
|
##args is now a list of character vectors
|
|
8
|
## First check to see if arguments are passed.
|
|
9
|
## if no args were given, print a warning and stop
|
|
10
|
if(length(args)==0) {stop("No parameters supplied, you must pass parameters")}
|
|
11
|
|
|
12
|
## Then cycle through each element of the list and evaluate the expressions.
|
|
13
|
eval(parse(text=args))
|
|
14
|
## now there is an i that corresponds to the row in the notdone object that will be processed.
|
|
15
|
|
|
16
|
## load libraries
|
|
17
|
require(sp)
|
|
18
|
require(rgdal)
|
|
19
|
require(reshape)
|
|
20
|
require(ncdf4)
|
|
21
|
require(geosphere)
|
|
22
|
require(raster)
|
|
23
|
require(spgrass6)
|
|
24
|
|
|
25
|
## specify some working directories
|
|
26
|
setwd("/nobackupp1/awilso10/mod06")
|
|
27
|
tempdir=tempdir() # to hold intermediate files
|
|
28
|
outdir="2_daily" # final daily product
|
|
29
|
tile="h11v08" #can move this to submit script if needed
|
|
30
|
## read in list of all files
|
|
31
|
load("allfiles.Rdata")
|
|
32
|
load("notdone.Rdata")
|
|
33
|
load("vars.Rdata")
|
|
34
|
|
|
35
|
date=notdone[i]
|
|
36
|
#file=fs$path[fs$dateid==date][1]
|
|
37
|
|
|
38
|
## print out some info
|
|
39
|
print(paste("Processing date ",date," for tile",tile))
|
|
40
|
|
|
41
|
## identify tile of interest
|
|
42
|
tile_bb=tb[tb$tile==tile,]
|
|
43
|
|
|
44
|
#### Function to generate hegtool parameter file for multi-band HDF-EOS file
|
|
45
|
swath2grid=function(file,vars,upleft,lowright){
|
|
46
|
print(paste("Starting file",basename(file)))
|
|
47
|
outfile=paste(tempdir(),"/",basename(file),sep="")
|
|
48
|
### First write the parameter file (careful, heg is very finicky!)
|
|
49
|
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
50
|
grp=paste("
|
|
51
|
BEGIN
|
|
52
|
INPUT_FILENAME=",file,"
|
|
53
|
OBJECT_NAME=mod06
|
|
54
|
FIELD_NAME=",vars$variable,"|
|
|
55
|
BAND_NUMBER = 1
|
|
56
|
OUTPUT_PIXEL_SIZE_X=1000
|
|
57
|
OUTPUT_PIXEL_SIZE_Y=1000
|
|
58
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
|
59
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
|
60
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
|
61
|
RESAMPLING_TYPE =NN
|
|
62
|
OUTPUT_PROJECTION_TYPE = SIN
|
|
63
|
#UTM_ZONE = 10
|
|
64
|
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 )
|
|
65
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
|
66
|
ELLIPSOID_CODE = WGS84
|
|
67
|
OUTPUT_TYPE = HDFEOS
|
|
68
|
OUTPUT_FILENAME = ",outfile,"
|
|
69
|
END
|
|
70
|
|
|
71
|
",sep="")
|
|
72
|
|
|
73
|
## if any remnants from previous runs remain, delete them
|
|
74
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
|
75
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
|
76
|
## write it to a file
|
|
77
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
|
78
|
## now run the swath2grid tool
|
|
79
|
## write the tiff file
|
|
80
|
log=system(paste("/nobackupp4/pvotava/software/heg/bin/swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
|
81
|
## clean up temporary files in working directory
|
|
82
|
# file.remove(paste("filetable.temp_",pid,sep=""))
|
|
83
|
print(log)
|
|
84
|
print(paste("Finished ", file))
|
|
85
|
}
|
|
86
|
|
|
87
|
#### Run the gridding procedure
|
|
88
|
lapply(fs$path[fs$dateid==date],swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max))
|
|
89
|
#upleft=paste(roi_bb[2,2],roi_bb[1,1]),lowright=paste(roi_bb[2,1],roi_bb[1,2]))
|
|
90
|
|
|
91
|
|
|
92
|
##############################################################
|
|
93
|
### Import to GRASS for processing
|
|
94
|
|
|
95
|
#fs$grass=paste(fs$month,fs$year,fs$file,sep="_")
|
|
96
|
#td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""))
|
|
97
|
|
|
98
|
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27/",pattern=paste(tile,".*hdf$",sep=""),full=T)
|
|
99
|
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
|
100
|
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 "
|
|
101
|
|
|
102
|
## fucntion to convert binary to decimal to assist in identifying correct values
|
|
103
|
b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
|
104
|
## for example:
|
|
105
|
b2d(c(T,T))
|
|
106
|
|
|
107
|
### create (or connect to) grass location
|
|
108
|
gisLocation="tmp"
|
|
109
|
gisMapset="mod06"
|
|
110
|
## set Grass to overwrite
|
|
111
|
Sys.setenv(GRASS_OVERWRITE=1)
|
|
112
|
Sys.setenv(DEBUG=1)
|
|
113
|
Sys.setenv(GRASS_GUI="txt")
|
|
114
|
|
|
115
|
## temporary objects to test function below
|
|
116
|
# i=1
|
|
117
|
#file=paste(outdir,"/",fs$file[1],sep="")
|
|
118
|
#date=as.Date("2000-05-23")
|
|
119
|
|
|
120
|
#.GRASS_CACHE=spgrass6:::.GRASS_CACHE
|
|
121
|
|
|
122
|
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
|
123
|
loadcloud<-function(date,fs){
|
|
124
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")
|
|
125
|
print(paste("Set up temporary grass session in",tf))
|
|
126
|
|
|
127
|
## set up temporary grass instance for this PID
|
|
128
|
initGRASS(gisBase="/nobackupp1/awilso10/software/grass-6.4.3svn",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
129
|
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""))
|
|
130
|
|
|
131
|
## Define region by importing one MOD11A1 raster.
|
|
132
|
print("Import one MOD11A1 raster to define grid")
|
|
133
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
134
|
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
135
|
system("g.region rast=modisgrid save=roi --overwrite")
|
|
136
|
|
|
137
|
## Identify which files to process
|
|
138
|
tfs=fs$file[fs$dateid==date]
|
|
139
|
nfs=length(tfs)
|
|
140
|
|
|
141
|
## loop through scenes and process QA flags
|
|
142
|
for(i in 1:nfs){
|
|
143
|
file=paste(tempdir,"/",tfs[i],sep="")
|
|
144
|
## Cloud Mask
|
|
145
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
146
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
147
|
## extract cloudy and 'confidently clear' pixels
|
|
148
|
system(paste("r.mapcalc <<EOF
|
|
149
|
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0
|
|
150
|
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3
|
|
151
|
EOF",sep=""))
|
|
152
|
|
|
153
|
## QA
|
|
154
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
|
155
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
156
|
## QA_CER
|
|
157
|
system(paste("r.mapcalc <<EOF
|
|
158
|
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
|
159
|
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3
|
|
160
|
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0
|
|
161
|
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1
|
|
162
|
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3
|
|
163
|
EOF",sep=""))
|
|
164
|
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1
|
|
165
|
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3
|
|
166
|
|
|
167
|
## Optical Thickness
|
|
168
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
|
169
|
output=paste("COT_",i,sep=""),
|
|
170
|
title="cloud_effective_radius",
|
|
171
|
flags=c("overwrite","o")) ; print("")
|
|
172
|
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
|
173
|
## keep only positive COT values where quality is 'useful' and 'very good' & scale to real units
|
|
174
|
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=""))
|
|
175
|
## set COT to 0 in clear-sky pixels
|
|
176
|
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
|
177
|
|
|
178
|
## Effective radius ##
|
|
179
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
|
180
|
output=paste("CER_",i,sep=""),
|
|
181
|
title="cloud_effective_radius",
|
|
182
|
flags=c("overwrite","o")) ; print("")
|
|
183
|
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
|
184
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
|
185
|
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
186
|
## set CER to 0 in clear-sky pixels
|
|
187
|
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
|
188
|
|
|
189
|
## Cloud Water Path
|
|
190
|
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
|
191
|
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
|
192
|
# flags=c("overwrite","o")) ; print("")
|
|
193
|
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
|
194
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
|
195
|
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
196
|
## set CER to 0 in clear-sky pixels
|
|
197
|
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
|
198
|
|
|
199
|
} #end loop through sub daily files
|
|
200
|
|
|
201
|
#### Now generate daily averages (or maximum in case of cloud flag)
|
|
202
|
|
|
203
|
system(paste("r.mapcalc <<EOF
|
|
204
|
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
205
|
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
206
|
COT_daily=COT_numer/COT_denom
|
|
207
|
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
208
|
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
209
|
CER_daily=CER_numer/CER_denom
|
|
210
|
CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),")
|
|
211
|
EOF",sep=""))
|
|
212
|
|
|
213
|
#### Write the files to a geotiff
|
|
214
|
execGRASS("r.out.gdal",input="CER_daily",output=paste(outdir,"/CER_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
215
|
execGRASS("r.out.gdal",input="COT_daily",output=paste(outdir,"/COT_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
216
|
execGRASS("r.out.gdal",input="CLD_daily",output=paste(outdir,"/CLD_",date,".tif",sep=""),nodata=99,flags=c("quiet"))
|
|
217
|
|
|
218
|
### delete the temporary files
|
|
219
|
unlink_.gislock()
|
|
220
|
system("/nobackupp1/awilso10/software/grass-6.4.3svn/etc/clean_temp")
|
|
221
|
system(paste("rm -R ",tf,sep=""))
|
|
222
|
### print update
|
|
223
|
}
|
|
224
|
|
|
225
|
|
|
226
|
###########################################
|
|
227
|
### Now run it
|
|
228
|
|
|
229
|
lapply(date,function(date) loadcloud(date,fs=fs))
|
|
230
|
|
|
231
|
|
|
232
|
|
|
233
|
print(paste(" ################################################################### Finished ",date,"
|
|
234
|
################################################################"))
|
|
235
|
|
|
236
|
|
|
237
|
## quit R
|
|
238
|
q("no")
|
|
239
|
|
|
240
|
|
