Revision 95354b03
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 196 | 196 |
plot=F |
| 197 | 197 |
if(plot){
|
| 198 | 198 |
i=1 |
| 199 |
system(paste("gdalinfo ",swaths[1]))
|
|
| 200 |
d=brick(lapply(outfiles,function(r) raster(r)))
|
|
| 201 |
plot(d)
|
|
| 199 |
system(paste("gdalinfo ",outfiles[19]))
|
|
| 200 |
d=lapply(outfiles,function(r) raster(r))
|
|
| 201 |
summary(d[[5]])
|
|
| 202 | 202 |
} |
| 203 | 203 |
#system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
| 204 | 204 |
|
| ... | ... | |
| 262 | 262 |
## Sensor Zenith ## extract first bit to keep only "low angle" observations |
| 263 | 263 |
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
| 264 | 264 |
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 265 |
|
|
| 266 |
## check for interpolation artefacts |
|
| 267 |
# execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))
|
|
| 268 |
# execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))
|
|
| 269 |
# execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))
|
|
| 270 |
|
|
| 271 |
## write out the table of weights for use in the neighborhood analysis to identify bad pixels from swtif |
|
| 272 |
p=75 #must be odd |
|
| 273 |
mat=matrix(rep(0,p*p),nrow=p) |
|
| 274 |
mat[0.5+p/2,]=1 |
|
| 275 |
cat(mat,file="weights.txt") |
|
| 276 |
execGRASS("r.neighbors",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"clump",sep=""),method="range",size=p,weight="weights.txt") # too slow!
|
|
| 277 |
system(paste("r.mapcalc \"SZ_",i,"_clump2=SZ_",i,"clump==0\"",sep=""))
|
|
| 278 |
|
|
| 279 |
# p=-50:50 |
|
| 280 |
# 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=""))
|
|
| 281 |
# 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=""))
|
|
| 282 |
# vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))
|
|
| 283 |
# colnames(vals)=c("value","count")
|
|
| 284 |
# vals$count=as.numeric(as.character(vals$count)) |
|
| 285 |
# vals$value=as.numeric(as.character(vals$value)) |
|
| 286 |
# vals=na.omit(vals) |
|
| 287 |
# vals$count[vals$value==1&vals$count>10] |
|
| 288 |
# |
|
| 289 |
#plot(p~value,data=vals) |
|
| 290 |
# print(sum(vals$p[vals$p>.1])) |
|
| 291 |
|
|
| 292 | 265 |
## Solar Zenith ## extract first bit to keep only "low angle" observations |
| 293 | 266 |
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
| 294 | 267 |
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| ... | ... | |
| 296 | 269 |
system(paste("r.mapcalc <<EOF
|
| 297 | 270 |
CM_fill_",i," = if(isnull(CM1_",i,"),1,0) |
| 298 | 271 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0) |
| 299 |
SZ_low_",i," = if(SZ_",i,"_clump2==0&SZ_",i,"<6000,1,0)
|
|
| 272 |
SZ_low_",i," = if(SZ_",i,"<6000,1,0) |
|
| 300 | 273 |
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0) |
| 301 | 274 |
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0) |
| 302 | 275 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null()) |
| 303 |
SZday_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==1,SZ_",i,",null())
|
|
| 304 |
SZnight_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==0,SZ_",i,",null())
|
|
| 276 |
SZday_",i," = if(CM_dayflag_",i,"==1,SZ_",i,",null()) |
|
| 277 |
SZnight_",i," = if(CM_dayflag_",i,"==0,SZ_",i,",null()) |
|
| 305 | 278 |
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null()) |
| 306 | 279 |
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null()) |
| 307 | 280 |
EOF",sep="")) |
| ... | ... | |
| 315 | 288 |
d2=stack( |
| 316 | 289 |
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
| 317 | 290 |
raster(readRAST6(paste("CM1_",i,sep=""))),
|
| 318 |
# raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
| 319 |
# raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
| 320 |
# raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
| 321 |
# raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
| 291 |
raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
| 292 |
raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
| 293 |
raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
| 294 |
raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
| 322 | 295 |
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
| 323 | 296 |
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
| 324 |
raster(readRAST6(paste("SZ_",i,sep=""))),
|
|
| 325 |
raster(readRAST6(paste("SZ_",i,"_clump",sep=""))),
|
|
| 326 |
raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))
|
|
| 297 |
raster(readRAST6(paste("SZ_",i,sep="")))
|
|
| 327 | 298 |
) |
| 328 | 299 |
plot(d2,add=F) |
| 329 | 300 |
} |
| ... | ... | |
| 346 | 317 |
if(plot){
|
| 347 | 318 |
ps=1:nfs |
| 348 | 319 |
ps=c(12,14,17) |
| 349 |
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))
|
|
| 350 |
sz_clump=brick(lapply(ps,function(i) raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))))
|
|
| 351 |
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
| 320 |
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZday_",i,sep="")))))
|
|
| 321 |
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMday_",i,sep="")))))
|
|
| 352 | 322 |
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
| 353 | 323 |
library(rasterVis) |
| 354 | 324 |
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100)) |
| 355 |
levelplot(sz_clump) |
|
| 356 | 325 |
levelplot(d) |
| 357 | 326 |
levelplot(d2) |
| 358 | 327 |
} |
| climate/procedures/Pleiades_MOD35.R | ||
|---|---|---|
| 1 | 1 |
#### Script to facilitate processing of MOD06 data |
| 2 |
|
|
| 2 |
### This script is meant to be run iteratively, rather than unsupervised. There are several steps that require manual checking (such as choosing the number of cores, etc.) |
|
| 3 |
|
|
| 4 |
## working directory |
|
| 3 | 5 |
setwd("/nobackupp1/awilso10/mod35")
|
| 4 | 6 |
|
| 7 |
## load libraries |
|
| 5 | 8 |
library(rgdal) |
| 6 | 9 |
library(raster) |
| 7 | 10 |
library(RSQLite) |
| 8 | 11 |
|
| 9 |
|
|
| 12 |
## flag to increase verbosity of output |
|
| 10 | 13 |
verbose=T |
| 11 | 14 |
|
| 12 | 15 |
## get MODLAND tile information |
| ... | ... | |
| 19 | 22 |
## Choose some tiles to process |
| 20 | 23 |
### list of tiles to process |
| 21 | 24 |
tiles=c("h10v08","h11v08","h12v08","h10v07","h11v07","h12v07") # South America
|
| 22 |
## a northern block of tiles |
|
| 25 |
## or a northern block of tiles
|
|
| 23 | 26 |
tiles=apply(expand.grid(paste("h",11:17,sep=""),v=c("v00","v01","v02","v03","v04")),1,function(x) paste(x,collapse="",sep=""))
|
| 24 | 27 |
## subset to MODLAND tiles |
| 25 | 28 |
alltiles=system("ls -r MODTILES/ | grep tif$ | cut -c1-6 | sort | uniq - ",intern=T)
|
| ... | ... | |
| 99 | 102 |
|
| 100 | 103 |
### report on what has already been processed |
| 101 | 104 |
print(paste(sum(!proclist$done)," out of ",nrow(proclist)," (",round(100*sum(!proclist$done)/nrow(proclist),2),"%) remain"))
|
| 105 |
stem(table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done])) |
|
| 102 | 106 |
table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done]) |
| 103 | 107 |
table(table(tile=proclist$tile[!proclist$done],year=proclist$year[!proclist$done])) |
| 104 | 108 |
|
| ... | ... | |
| 112 | 116 |
tp=T # rerun everything |
| 113 | 117 |
tp=((!proclist$done)&proclist$avail) #date-tiles to process |
| 114 | 118 |
table(Available=proclist$avail,Completed=proclist$done) |
| 119 |
table(tp) |
|
| 115 | 120 |
|
| 116 | 121 |
write.table(paste("--verbose ",script," --date ",proclist$date[tp]," --verbose T --tile ",proclist$tile[tp],sep=""),
|
| 117 | 122 |
file=paste("notdone.txt",sep=""),row.names=F,col.names=F,quote=F)
|
| 118 | 123 |
|
| 124 |
## try running it once for a single tile-date to get estimate of time/tile-day |
|
| 125 |
test=F |
|
| 126 |
if(test){
|
|
| 127 |
i=2 |
|
| 128 |
time1=system.time(system(paste("Rscript --verbose ",script," --date ",proclist$date[i]," --verbose T --tile ",proclist$tile[i],sep="")))
|
|
| 129 |
hours=round(length(proclist$date[tp])*142/60/60) |
|
| 130 |
hours=round(length(proclist$date[tp])*time1[3]/60/60,1) |
|
| 131 |
hours/240 |
|
| 132 |
print(paste("Based on runtime of previous command, it will take",hours," hours to process the full set"))
|
|
| 133 |
} |
|
| 134 |
|
|
| 135 |
|
|
| 119 | 136 |
### qsub script |
| 120 | 137 |
cat(paste("
|
| 121 | 138 |
#PBS -S /bin/bash |
| 122 |
#PBS -l select=20:ncpus=8:mpiprocs=8
|
|
| 139 |
#PBS -l select=28:ncpus=8:mpiprocs=8
|
|
| 123 | 140 |
##PBS -l select=100:ncpus=8:mpiprocs=8 |
| 124 | 141 |
##PBS -l walltime=8:00:00 |
| 125 |
#PBS -l walltime=4:00:00
|
|
| 142 |
#PBS -l walltime=2:00:00
|
|
| 126 | 143 |
#PBS -j n |
| 127 | 144 |
#PBS -m be |
| 128 | 145 |
#PBS -N mod35 |
| 129 |
#PBS -q normal |
|
| 130 |
##PBS -q devel
|
|
| 146 |
##PBS -q normal
|
|
| 147 |
#PBS -q devel |
|
| 131 | 148 |
#PBS -V |
| 132 | 149 |
|
| 133 | 150 |
#CORES=800 |
| 134 |
CORES=160
|
|
| 151 |
CORES=224
|
|
| 135 | 152 |
|
| 136 | 153 |
HDIR=/u/armichae/pr/ |
| 137 | 154 |
source $HDIR/etc/environ.sh |
| ... | ... | |
| 152 | 169 |
system(paste("cat notdone.txt | head",sep=""))
|
| 153 | 170 |
system(paste("cat notdone.txt | wc -l ",sep=""))
|
| 154 | 171 |
|
| 172 |
|
|
| 155 | 173 |
## Submit it |
| 156 | 174 |
system(paste("qsub mod35_qsub",sep=""))
|
| 157 | 175 |
|
Also available in: Unified diff
Updated to use new swtif and remove additional QC that identified strips introduced by old swtif. Some tiling artifacts still present.