Revision 848e34fd
Added by Adam Wilson about 11 years ago
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 24 | 24 |
'date', 'd', 1, 'character', |
| 25 | 25 |
'tile', 't', 1, 'character', |
| 26 | 26 |
'verbose','v',1,'logical', |
| 27 |
'profile','p',0,'logical', |
|
| 27 | 28 |
'help', 'h', 0, 'logical' |
| 28 | 29 |
), ncol=4, byrow=TRUE)) |
| 29 | 30 |
if ( !is.null(opta$help) ) |
| ... | ... | |
| 36 | 37 |
testing=F |
| 37 | 38 |
platform="pleiades" |
| 38 | 39 |
|
| 40 |
## record profiling information if requested |
|
| 41 |
if(opta$profile) Rprof("/nobackupp1/awilso10/mod35/log/profile.out")
|
|
| 42 |
|
|
| 39 | 43 |
## default date and tile to play with (will be overwritten below when running in batch) |
| 40 | 44 |
if(testing){
|
| 41 | 45 |
date="20090129" |
| ... | ... | |
| 95 | 99 |
|
| 96 | 100 |
|
| 97 | 101 |
### print some status messages |
| 98 |
if(verbose) print(paste("Processing tile",tile," for date",date))
|
|
| 102 |
if(verbose) writeLines(paste("STATUS: Beginning ",tile,date))
|
|
| 99 | 103 |
|
| 100 | 104 |
## load tile information and get bounding box |
| 101 | 105 |
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata") |
| ... | ... | |
| 128 | 132 |
## find the swaths on disk (using datadir) |
| 129 | 133 |
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T) |
| 130 | 134 |
|
| 131 |
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk")) |
|
| 135 |
### print some status messages |
|
| 136 |
if(verbose) writeLines(paste("STATUS:swaths tile:",tile,"date:",date,"swathIDs:",nrow(fs)," swathsOnDisk:",length(swaths)))
|
|
| 132 | 137 |
|
| 133 | 138 |
|
| 134 | 139 |
## define function that grids swaths |
| ... | ... | |
| 185 | 190 |
|
| 186 | 191 |
|
| 187 | 192 |
############################################################################# |
| 188 |
## Check output |
|
| 189 |
## confirm at least one file for this date is present. If not, quit. |
|
| 190 |
outfiles=list.files(tempdir(),full=T,pattern=paste(basename(swaths),"$",sep="",collapse="|")) |
|
| 191 |
if(!any(file.exists(outfiles))) {
|
|
| 193 |
## check for zero dimension in HDFs |
|
| 194 |
## occasionlly swtif will output a hdf with a resolution of 0. Not sure why, but drop them here. |
|
| 195 |
CMcheck=list.files(pattern="CM_.*hdf$") #list of files to check |
|
| 196 |
CM_0=do.call(c,lapply(CMcheck, function(f) any(res(raster(f))==0))) |
|
| 197 |
keep=sub("CM_","",CMcheck[!CM_0])
|
|
| 198 |
if(length(keep)<length(CMcheck)){writeLines(paste("Warning (Resolution of zero): ",paste(sub("CM_","",CMcheck)[!sub("CM_","",CMcheck)%in%keep],collapse=",")," from ",tile," for ",date))}
|
|
| 199 |
outfiles=list.files(tempdir(),full=T,pattern=paste(keep,"$",sep="",collapse="|")) |
|
| 200 |
if(length(outfiles)==0) {
|
|
| 192 | 201 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
| 193 | 202 |
q("no",status=0)
|
| 194 | 203 |
} |
| 195 | 204 |
|
| 205 |
## confirm at least one file for this date is present. If not, quit. |
|
| 206 |
#outfiles=list.files(tempdir(),full=T,pattern=paste(basename(swaths),"$",sep="",collapse="|")) |
|
| 207 |
#if(!any(file.exists(outfiles))) {
|
|
| 208 |
# print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
|
| 209 |
# q("no",status=0)
|
|
| 210 |
#} |
|
| 211 |
|
|
| 212 |
|
|
| 196 | 213 |
plot=F |
| 197 | 214 |
if(plot){
|
| 198 | 215 |
i=1 |
| ... | ... | |
| 251 | 268 |
|
| 252 | 269 |
## loop through scenes and process QA flags |
| 253 | 270 |
for(i in 1:nfs){
|
| 254 |
bfile=tfs[i]
|
|
| 271 |
bfile=tfs[i] |
|
| 255 | 272 |
## Read in the data from the HDFs |
| 256 | 273 |
## Cloud Mask |
| 274 |
GDALinfo(paste("CM_",bfile,sep=""),returnStats=F,silent=T)
|
|
| 257 | 275 |
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
| 258 | 276 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 259 | 277 |
## QA ## extract first bit to keep only "useful" values of cloud mask |
| ... | ... | |
| 399 | 417 |
# system(paste("rm -frR ",tempdir(),sep=""))
|
| 400 | 418 |
|
| 401 | 419 |
|
| 402 |
## print out some info |
|
| 403 |
print(paste("####################### Finished ",tile,"-",date, "###################################"))
|
|
| 420 |
### print some status messages |
|
| 421 |
if(verbose) writeLines(paste("STATUS:end tile:",tile,"date:",date,"swathIDs:",nrow(fs)," swathsOnDisk:",length(swaths),"fileExists:",file.exists(ncfile)))
|
|
| 422 |
|
|
| 423 |
## turn off the profiler |
|
| 424 |
if(opta$profile) Rprof(NULL) |
|
| 404 | 425 |
|
| 405 |
## delete old files |
|
| 406 |
#system("cleartemp")
|
|
| 407 | 426 |
|
| 408 | 427 |
## quit |
| 409 | 428 |
q("no",status=0)
|
| climate/procedures/Pleiades_MOD35.R | ||
|---|---|---|
| 27 | 27 |
## subset to MODLAND tiles |
| 28 | 28 |
alltiles=system("ls -r MODTILES/ | grep tif$ | cut -c1-6 | sort | uniq - ",intern=T)
|
| 29 | 29 |
|
| 30 |
## or run all tiles |
|
| 31 |
#tiles=alltiles |
|
| 32 |
|
|
| 30 | 33 |
## subset to tiles in global region (not outside global boundary in sinusoidal projection) |
| 31 | 34 |
tiles=tiles[tiles%in%alltiles] |
| 32 | 35 |
|
| ... | ... | |
| 38 | 41 |
|
| 39 | 42 |
outdir="daily/" #paste("daily/",tile,sep="")
|
| 40 | 43 |
|
| 41 |
##find swaths in region from sqlite database for the specified date/tile
|
|
| 44 |
##find swaths in region from sqlite database for the specified tile |
|
| 42 | 45 |
## this takes a while, about 30 minutes, so only rebuild if you need to update what's available... |
| 43 | 46 |
rebuildswathtable=F |
| 44 | 47 |
if(rebuildswathtable){
|
| ... | ... | |
| 69 | 72 |
save(fs,swaths,file="swathtile.Rdata") |
| 70 | 73 |
} |
| 71 | 74 |
|
| 72 |
load("swathtile.Rdata")
|
|
| 75 |
if(!exists("fs")) load("swathtile.Rdata")
|
|
| 73 | 76 |
|
| 74 | 77 |
if(verbose) print(paste("###############",nrow(fs)," swath IDs recieved from database"))
|
| 75 | 78 |
|
| ... | ... | |
| 118 | 121 |
table(Available=proclist$avail,Completed=proclist$done) |
| 119 | 122 |
table(tp) |
| 120 | 123 |
|
| 121 |
write.table(paste("--verbose ",script," --date ",proclist$date[tp]," --verbose T --tile ",proclist$tile[tp],sep=""),
|
|
| 124 |
write.table(paste("--verbose ",script," --date ",proclist$date[tp]," --verbose T --profile F --tile ",proclist$tile[tp],sep=""),
|
|
| 122 | 125 |
file=paste("notdone.txt",sep=""),row.names=F,col.names=F,quote=F)
|
| 123 | 126 |
|
| 124 | 127 |
## try running it once for a single tile-date to get estimate of time/tile-day |
| 125 | 128 |
test=F |
| 126 | 129 |
if(test){
|
| 127 | 130 |
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) |
|
| 131 |
time1=system.time(system(paste("Rscript --verbose ",script," --date ",proclist$date[i]," --profile T --verbose T --tile ",proclist$tile[i],sep="")))
|
|
| 132 |
hours=round(length(proclist$date[tp])*142/60/60); hours
|
|
| 130 | 133 |
hours=round(length(proclist$date[tp])*time1[3]/60/60,1); hours |
| 131 |
hours/400 |
|
| 132 |
print(paste("Based on runtime of previous command, it will take",hours," hours to process the full set"))
|
|
| 134 |
nodes=100 |
|
| 135 |
threads=nodes*8 |
|
| 136 |
writeLines(paste(" ################### \n Hours per date-tile:",round(time1[3]/60/60,2),"\n Date-tiles to process:",sum(tp)," \n Estimated CPU time: ",hours,"hours \n With ",threads,"threads:",round(hours/threads,2),"hours \n ###################"))
|
|
| 133 | 137 |
} |
| 134 | 138 |
|
| 139 |
### Set up submission script |
|
| 140 |
queue="devel" |
|
| 141 |
queue="normal" #"devel" |
|
| 142 |
nodes=50 |
|
| 143 |
walltime=2 |
|
| 135 | 144 |
|
| 136 |
### qsub script
|
|
| 145 |
### write qsub script to disk
|
|
| 137 | 146 |
cat(paste("
|
| 138 | 147 |
#PBS -S /bin/bash |
| 139 |
#PBS -l select=50:ncpus=8:mpiprocs=8 |
|
| 140 |
##PBS -l select=100:ncpus=8:mpiprocs=8 |
|
| 141 |
##PBS -l walltime=8:00:00 |
|
| 142 |
#PBS -l walltime=2:00:00 |
|
| 148 |
#PBS -l select=",nodes,":ncpus=8:mpiprocs=8 |
|
| 149 |
#PBS -l walltime=",walltime,":00:00 |
|
| 143 | 150 |
#PBS -j n |
| 144 | 151 |
#PBS -m be |
| 145 | 152 |
#PBS -N mod35 |
| 146 |
##PBS -q normal |
|
| 147 |
#PBS -q devel |
|
| 153 |
#PBS -q ",queue," |
|
| 148 | 154 |
#PBS -V |
| 149 | 155 |
|
| 150 |
#CORES=800 |
|
| 151 |
CORES=400 |
|
| 152 |
|
|
| 156 |
CORES=",nodes*8," |
|
| 153 | 157 |
HDIR=/u/armichae/pr/ |
| 154 | 158 |
source $HDIR/etc/environ.sh |
| 155 | 159 |
source /u/awilso10/environ.sh |
| 156 | 160 |
source /u/awilso10/.bashrc |
| 157 | 161 |
IDIR=/nobackupp1/awilso10/mod35/ |
| 158 |
##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt |
|
| 159 | 162 |
WORKLIST=$IDIR/notdone.txt |
| 160 | 163 |
EXE=Rscript |
| 161 | 164 |
LOGSTDOUT=$IDIR/log/mod35_stdout |
| ... | ... | |
| 166 | 169 |
|
| 167 | 170 |
### Check the files |
| 168 | 171 |
system(paste("cat mod35_qsub",sep=""))
|
| 169 |
system(paste("cat notdone.txt | head",sep=""))
|
|
| 172 |
system(paste("cat notdone.txt | head -n 4",sep=""))
|
|
| 170 | 173 |
system(paste("cat notdone.txt | wc -l ",sep=""))
|
| 171 | 174 |
|
| 175 |
## start interactive job on compute node for debugging |
|
| 176 |
# system("qsub -I -l walltime=2:00:00 -lselect=2:ncpus=16:model=san -q devel")
|
|
| 177 |
|
|
| 172 | 178 |
|
| 173 | 179 |
## Submit it |
| 174 | 180 |
system(paste("qsub mod35_qsub",sep=""))
|
| ... | ... | |
| 204 | 210 |
file=paste("notdone_climate.txt",sep=""),row.names=F,col.names=F,quote=F)
|
| 205 | 211 |
|
| 206 | 212 |
## delay start until previous jobs have finished? |
| 207 |
delay=T
|
|
| 213 |
delay=F
|
|
| 208 | 214 |
## check running jobs to get JobID of job you want to wait for |
| 209 | 215 |
system("qstat -u awilso10",intern=T)
|
| 210 | 216 |
## enter JobID here: |
| ... | ... | |
| 250 | 256 |
## check progress |
| 251 | 257 |
system("qstat -u awilso10")
|
| 252 | 258 |
|
| 253 |
## start interactive job on compute node for debugging |
|
| 254 |
# system("qsub -I -l walltime=2:00:00 -lselect=2:ncpus=16:model=san -q devel")
|
|
| 255 | 259 |
|
| 256 | 260 |
|
| 257 | 261 |
################################################################# |
| 258 | 262 |
### copy the files back to Yale |
| 259 | 263 |
|
| 260 | 264 |
|
| 261 |
system("ssh lou")
|
|
| 265 |
#system("ssh lou")
|
|
| 262 | 266 |
#scp `find MOD35/summary -name "MOD35_h[0-9][0-9]v[0-9][0-9].nc"` adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/summary/ |
| 263 | 267 |
system("rsync -cavv `find summary -name \"MOD35_h[0-9][0-9]v[0-9][0-9]_mean.nc\"` adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/summary/")
|
| 264 | 268 |
system("rsync -cavv `find summary -name \"MOD35_h[0-9][0-9]v[0-9][0-9].nc\"` adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/summary/")
|
Also available in: Unified diff
Submitted global run for 2009