| 40 |
40 |
gdir="output/"
|
| 41 |
41 |
datadir="data/modis/MOD06_L2_hdf"
|
| 42 |
42 |
outdir="data/modis/MOD06_L2_hdf2"
|
| 43 |
|
|
|
43 |
tifdir="data/modis/MOD06_L2_tif"
|
|
44 |
|
| 44 |
45 |
fs=data.frame(
|
| 45 |
46 |
path=list.files(datadir,full=T,recursive=T,pattern="hdf"),
|
| 46 |
47 |
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))
|
| ... | ... | |
| 123 |
124 |
paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
| 124 |
125 |
print(paste("Finished ", file))
|
| 125 |
126 |
}
|
| 126 |
|
|
|
127 |
|
| 127 |
128 |
|
| 128 |
129 |
### update fs with completed files
|
| 129 |
130 |
fs$complete=fs$file%in%list.files(outdir,pattern="hdf$")
|
| ... | ... | |
| 156 |
157 |
gisDbase="/media/data/grassdata"
|
| 157 |
158 |
gisLocation="oregon"
|
| 158 |
159 |
gisMapset="mod06"
|
|
160 |
## set Grass to overwrite
|
|
161 |
Sys.setenv(GRASS_OVERWRITE=1)
|
|
162 |
Sys.setenv(DEBUG=0)
|
| 159 |
163 |
|
| 160 |
|
initGRASS(gisBase="/usr/lib/grass64",gisDbase=gisDbase,SG=td,override=T,
|
| 161 |
|
location=gisLocation,mapset=gisMapset)
|
|
164 |
initGRASS(gisBase="/usr/lib/grass64",SG=td,gisDbase=gisDbase,location=gisLocation,mapset="PERMANENT",override=T,pid=Sys.getpid())
|
| 162 |
165 |
getLocationProj()
|
|
166 |
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=""))
|
| 163 |
167 |
|
| 164 |
|
system("g.mapset mapset=PERMANENT")
|
|
168 |
#system("g.mapset PERMANENT")
|
| 165 |
169 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),
|
| 166 |
170 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
| 167 |
171 |
system("g.region rast=modisgrid save=roi --overwrite")
|
| 168 |
|
system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs\"",sep=""))
|
| 169 |
172 |
system("g.region roi")
|
| 170 |
173 |
system("g.region -p")
|
|
174 |
getLocationProj()
|
|
175 |
|
| 171 |
176 |
|
| 172 |
|
i=1
|
|
177 |
i=1
|
| 173 |
178 |
file=paste(outdir,"/",fs$file[1],sep="")
|
| 174 |
179 |
date=as.Date("2000-03-02")
|
| 175 |
180 |
|
| 176 |
181 |
loadcloud<-function(date,fs){
|
| 177 |
|
sink(file=paste("output/",format(date,"%Y%m%d"),".txt",sep=""),type="output")
|
|
182 |
|
| 178 |
183 |
## Identify which files to process
|
| 179 |
184 |
tfs=fs$file[fs$date==date]
|
| 180 |
185 |
nfs=length(tfs)
|
| 181 |
|
|
| 182 |
|
## create new mapset to hold all data for this day
|
| 183 |
|
system(paste("g.mapset --quiet -c mapset=",gisMapset,"_",format(date,"%Y%m%d"),sep=""))
|
| 184 |
|
# file.copy(paste(gisDbase,"/",gisLocation,"/PERMANENT/DEFAULT_WIND",sep=""),paste(gisDbase,"/",gisLocation,"/",gisMapset,"_",format(date,"%Y%m%d"),"/WIND",sep=""))
|
|
186 |
unlink_.gislock()
|
|
187 |
## set new PID for this grass process (running within a spawned R session if using multicore)
|
|
188 |
set.GIS_LOCK(Sys.getpid())
|
|
189 |
## create new mapset to hold all data for this day
|
|
190 |
system(paste("g.mapset -c mapset=",gisMapset,"_",format(date,"%Y%m%d"),sep=""))
|
|
191 |
# file.copy(paste(gisDbase,"/",gisLocation,"/PERMANENT/DEFAULT_WIND",sep=""),paste(gisDbase,"/",gisLocation,"/",gisMapset,"_",format(date,"%Y%m%d"),"/WIND",sep=""))
|
| 185 |
192 |
system("g.region roi@PERMANENT")
|
| 186 |
|
|
|
193 |
print(date)
|
|
194 |
print(gmeta6())
|
|
195 |
print(Sys.getpid())
|
| 187 |
196 |
## loop through scenes and process QA flags
|
| 188 |
197 |
for(i in 1:nfs){
|
| 189 |
|
file=paste(outdir,"/",tfs[i],sep="")
|
| 190 |
|
## Cloud Mask
|
| 191 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
| 192 |
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o","quiet")) ; print("")
|
|
198 |
file=paste(outdir,"/",tfs[i],sep="")
|
|
199 |
## Cloud Mask
|
|
200 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
201 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 193 |
202 |
## extract cloudy and 'confidently clear' pixels
|
| 194 |
203 |
system(paste("r.mapcalc <<EOF
|
| 195 |
204 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0
|
| 196 |
205 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3
|
| 197 |
206 |
EOF",sep=""))
|
| 198 |
207 |
|
| 199 |
|
|
| 200 |
208 |
## QA
|
| 201 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| 202 |
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o","quiet")) ; print("")
|
|
209 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
|
210 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 203 |
211 |
## QA_CER
|
| 204 |
212 |
system(paste("r.mapcalc <<EOF
|
| 205 |
213 |
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
| ... | ... | |
| 244 |
252 |
} #end loop through sub daily files
|
| 245 |
253 |
|
| 246 |
254 |
#### Now generate daily averages
|
| 247 |
|
|
|
255 |
|
| 248 |
256 |
system(paste("r.mapcalc <<EOF
|
| 249 |
257 |
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
| 250 |
258 |
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
| ... | ... | |
| 255 |
263 |
EOF",sep=""))
|
| 256 |
264 |
|
| 257 |
265 |
#### Write the file to a geotiff
|
| 258 |
|
execGRASS("r.out.gdal",input="CER_daily",output=paste("data/modis/MOD06_L2_tif/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
| 259 |
|
execGRASS("r.out.gdal",input="COT_daily",output=paste("data/modis/MOD06_L2_tif/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
|
266 |
execGRASS("r.out.gdal",input="CER_daily",output=paste(tifdir,"/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
|
267 |
execGRASS("r.out.gdal",input="COT_daily",output=paste(tifdir,"/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
|
268 |
|
|
269 |
### delete the temporary files
|
|
270 |
unlink_.gislock()
|
|
271 |
system("/usr/lib/grass64/etc/clean_temp")
|
|
272 |
# system(paste("rm -R ",gmeta6()$GISDBASE,"/",gmeta6()$LOCATION_NAME,"/",gmeta6()$MAPSET,sep=""))
|
| 260 |
273 |
|
| 261 |
274 |
}
|
| 262 |
275 |
|
| ... | ... | |
| 264 |
277 |
###########################################
|
| 265 |
278 |
### Now run it
|
| 266 |
279 |
|
| 267 |
|
tdates=sort(unique(fs$date))
|
| 268 |
|
|
|
280 |
tdates=sort(unique(fs$date))
|
|
281 |
done=tdates%in%as.Date(substr(list.files("data/modis/MOD06_L2_tif"),5,12),"%Y%m%d")
|
|
282 |
table(done)
|
|
283 |
tdates=tdates[!done]
|
| 269 |
284 |
|
| 270 |
|
mclapply(tdates[1:10],loadcloud,fs=fs,mc.cores=1)
|
|
285 |
lapply(tdates,function(date) loadcloud(date,fs=fs))
|
| 271 |
286 |
|
| 272 |
287 |
|
| 273 |
288 |
|
| ... | ... | |
| 282 |
297 |
|
| 283 |
298 |
|
| 284 |
299 |
|
| 285 |
|
d=readRAST6("COT")
|
| 286 |
|
d$CER=readRAST6("CER")$CER
|
| 287 |
|
plot(COT~CER,data=d@data)
|
| 288 |
|
summary(lm(COT~CER,data=d@data))
|
| 289 |
|
|
|
300 |
#######################################################################################33
|
|
301 |
### Produce the monthly averages
|
| 290 |
302 |
|
|
303 |
## get list of daily files
|
|
304 |
fs2=data.frame(
|
|
305 |
path=list.files(outdir,full=T,recursive=T,pattern="hdf"),
|
|
306 |
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))
|
|
307 |
fs$date=as.Date(substr(fs$file,11,17),"%Y%j")
|
|
308 |
fs$month=format(fs$date,"%m")
|
|
309 |
fs$year=format(fs$date,"%Y")
|
|
310 |
fs$time=substr(fs$file,19,22)
|
|
311 |
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))
|
|
312 |
fs$dateid=format(fs$date,"%Y%m%d")
|
|
313 |
fs$path=as.character(fs$path)
|
|
314 |
fs$file=as.character(fs$file)
|
| 291 |
315 |
|
| 292 |
316 |
|
| 293 |
317 |
# read in data as single spatialgrid
|
| 294 |
318 |
ms=c("01","02","03","04","05","06","07","08","09","10","11","12")
|
| 295 |
|
for (m in ms){
|
|
319 |
|
|
320 |
mclapply(ms, function(m){
|
| 296 |
321 |
d=readRAST6(fs$grass[1])
|
| 297 |
322 |
projection(d)=projection(td)
|
| 298 |
323 |
d@data=as.data.frame(do.call(cbind,mclapply(which(fs$month==m),function(i){
|
First 10-year processing of daily mean tifs