Revision b3344197
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_Climatology.r | ||
|---|---|---|
| 10 | 10 |
'verbose','v',1,'logical' |
| 11 | 11 |
), ncol=4, byrow=TRUE)) |
| 12 | 12 |
|
| 13 |
tile=opta$tile #tile="h11v08"
|
|
| 13 |
tile=opta$tile #tile="h00v08"
|
|
| 14 | 14 |
verbose=opta$verbose #print out extensive information for debugging? |
| 15 | 15 |
|
| 16 | 16 |
## set working directory |
| 17 |
setwd("/u/awilso10/MOD35")
|
|
| 17 |
setwd("/nobackupp1/awilso10/mod35")
|
|
| 18 |
#setwd("/u/awilso10/MOD35")
|
|
| 18 | 19 |
|
| 19 | 20 |
### directory containing daily files |
| 20 | 21 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| ... | ... | |
| 68 | 69 |
#" -a missing_value,PClear,o,b,255 ", |
| 69 | 70 |
#" -a _FillValue,PClear,d,b,255 ", |
| 70 | 71 |
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ", |
| 71 |
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Daily Timeseries\" ",
|
|
| 72 |
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Summaries\" ",
|
|
| 72 | 73 |
" -a institution,global,o,c,\"Yale University\" ", |
| 73 |
" -a source,global,o,c,\"MODIS Cloud Mask (MOD35)\" ", |
|
| 74 |
" -a source,global,o,c,\"MODIS Collection 6 Cloud Mask (MOD35)\" ",
|
|
| 74 | 75 |
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ", |
| 75 | 76 |
outdir2,"/MOD35_",tile,"_daily.nc",sep="")) |
| 76 | 77 |
|
| ... | ... | |
| 91 | 92 |
|
| 92 | 93 |
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15) |
| 93 | 94 |
|
| 95 |
## Overall Means |
|
| 96 |
if(verbose) print(paste("Calculating the overall mean:",tile))
|
|
| 97 |
system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setmon,1 -setday,1 -mulc,-1 -subc,100 -timmean -selyear,2009 ",outdir2,"/MOD35_",tile,"_daily.nc ",outdir2,"/MOD35_",tile,"_2009mean.nc",sep=""),wait=T)
|
|
| 98 |
system(paste(ncopath,"ncrename -v PClear,PCloud ",outdir2,"/MOD35_",tile,"_2009mean.nc",sep="")) |
|
| 99 |
system(paste(ncopath,"ncatted ", |
|
| 100 |
" -a long_name,PCloud,o,c,\"Mean Probability of Cloud\" ", |
|
| 101 |
" -a missing_value,PCloud,o,b,255 ", |
|
| 102 |
" -a _FillValue,PCloud,d,b,255 ", |
|
| 103 |
outdir2,"/MOD35_",tile,"_2009mean.nc",sep="")) |
|
| 104 |
|
|
| 94 | 105 |
## Monthly means |
| 95 | 106 |
if(verbose) print(paste("Calculating the monthly means:",tile))
|
| 96 | 107 |
system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
| ... | ... | |
| 110 | 121 |
#system(paste("scp summary/MOD35_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))
|
| 111 | 122 |
#system(paste("ncdump -h ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",sep=""))
|
| 112 | 123 |
|
| 113 |
|
|
| 114 | 124 |
## Monthly standard deviation |
| 115 | 125 |
if(verbose) print(paste("Calculating the monthly SD:",tile))
|
| 116 | 126 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -mulc,-1 -subc,100 -monmean ",
|
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 35 | 35 |
|
| 36 | 36 |
|
| 37 | 37 |
## default date and tile to play with (will be overwritten below when running in batch) |
| 38 |
date="20000410" |
|
| 39 |
tile="h11v08" |
|
| 38 |
#date="20000410"
|
|
| 39 |
#tile="h11v08"
|
|
| 40 | 40 |
platform="pleiades" |
| 41 | 41 |
verbose=T |
| 42 | 42 |
|
| ... | ... | |
| 58 | 58 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
| 59 | 59 |
## specify working directory |
| 60 | 60 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 61 |
basedir="/nobackupp1/awilso10/mod35/tmp/" #directory to hold files temporarily before transferring to lou
|
|
| 61 |
basedir="/nobackupp1/awilso10/mod35/daily/" #directory to hold files temporarily before transferring to lou
|
|
| 62 | 62 |
setwd(tempdir()) |
| 63 | 63 |
## grass database |
| 64 | 64 |
gisBase="/u/armichae/pr/grass-6.4.2/" |
| ... | ... | |
| 207 | 207 |
## create output directory if needed |
| 208 | 208 |
## Identify output file |
| 209 | 209 |
ncfile=paste(basedir,"MOD35_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
| 210 |
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile)) |
|
| 210 |
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile,recursive=T))
|
|
| 211 | 211 |
|
| 212 | 212 |
## set up temporary grass instance for this PID |
| 213 | 213 |
if(verbose) print(paste("Set up temporary grass session in",tf))
|
| ... | ... | |
| 318 | 318 |
} |
| 319 | 319 |
|
| 320 | 320 |
############ copy files to lou |
| 321 |
if(platform=="pleiades"){
|
|
| 322 |
archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
|
| 323 |
system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
|
| 324 |
system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
|
| 325 |
file.remove(ncfile) |
|
| 326 |
file.remove(paste(ncfile,".aux.xml",sep="")) |
|
| 327 |
} |
|
| 321 |
#if(platform=="pleiades"){
|
|
| 322 |
# archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
|
| 323 |
# system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
|
| 324 |
# system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
|
| 325 |
# file.remove(ncfile)
|
|
| 326 |
# file.remove(paste(ncfile,".aux.xml",sep=""))
|
|
| 327 |
#}
|
|
| 328 | 328 |
|
| 329 | 329 |
|
| 330 | 330 |
### delete the temporary files |
| 331 |
unlink_.gislock() |
|
| 332 |
system(paste("rm -frR ",tempdir(),sep=""))
|
|
| 331 |
# unlink_.gislock()
|
|
| 332 |
# system(paste("rm -frR ",tempdir(),sep=""))
|
|
| 333 | 333 |
|
| 334 | 334 |
|
| 335 | 335 |
## print out some info |
| climate/procedures/Pleiades_MOD35.R | ||
|---|---|---|
| 19 | 19 |
## delete temporary log file that can grow to GB |
| 20 | 20 |
system("rm /nobackupp1/awilso10/software/heg/TOOLKIT_MTD/runtime/LogStatus")
|
| 21 | 21 |
|
| 22 |
|
|
| 23 |
tile="h11v08" # Venezuela |
|
| 24 |
#tile="h11v07" # Venezuela coast |
|
| 25 |
#tile="h09v04" # Oregon |
|
| 26 |
tile="h21v09" #Kenya |
|
| 27 |
|
|
| 28 | 22 |
### list of tiles to process |
| 29 |
tiles=c("h11v08","h21v09","h08v04","h09v04","h08v05","h09v05","h20v11","h31v11")
|
|
| 30 | 23 |
tiles=c("h10v08","h11v08","h12v08","h10v07","h11v07","h12v07") # South America
|
| 31 | 24 |
|
| 32 | 25 |
## subset to MODLAND tiles |
| 33 |
modlandtiles=system("ls -r /nobackupp4/datapool/modis/MOD11A1.005/2010* | grep hdf$ | cut -c18-23 | sort | uniq - ",intern=T)
|
|
| 26 |
modlandtiles=system("ls -r /nobackupp4/datapool/modis/MOD11A1.005/2010* | grep hdf$ | cut -c18-23 | sort | uniq - ",intern=T)
|
|
| 34 | 27 |
tb$land=tb$tile%in%modlandtiles |
| 35 | 28 |
tiles=tb$tile[tb$land] |
| 36 | 29 |
|
| ... | ... | |
| 43 | 36 |
outdir="daily/" #paste("daily/",tile,sep="")
|
| 44 | 37 |
|
| 45 | 38 |
##find swaths in region from sqlite database for the specified date/tile |
| 46 |
## path to swath database |
|
| 47 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 48 |
con=dbConnect("SQLite", dbname = db)
|
|
| 49 |
fs=do.call(rbind.data.frame,lapply(1:nrow(tile_bb),function(i){
|
|
| 50 |
d=dbGetQuery(con,paste("SELECT * from swath_geo6
|
|
| 39 |
## this takes a while, about 30 minutes, so only rebuild if you need to update what's available... |
|
| 40 |
rebuildswathtable=F |
|
| 41 |
if(rebuildswathtable){
|
|
| 42 |
## path to swath database |
|
| 43 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 44 |
con=dbConnect("SQLite", dbname = db)
|
|
| 45 |
fs=do.call(rbind.data.frame,lapply(1:nrow(tile_bb),function(i){
|
|
| 46 |
d=dbGetQuery(con,paste("SELECT * from swath_geo6
|
|
| 51 | 47 |
WHERE east>=",tile_bb$lon_min[i]," AND |
| 52 | 48 |
west<=",tile_bb$lon_max[i]," AND |
| 53 | 49 |
north>=",tile_bb$lat_min[i]," AND |
| 54 | 50 |
south<=",tile_bb$lat_max[i]) |
| 55 |
) |
|
| 56 |
d$tile=tile_bb$tile[i] |
|
| 57 |
print(paste("Finished tile",tile_bb$tile[i]))
|
|
| 58 |
return(d) |
|
| 59 |
})) |
|
| 60 |
con=dbDisconnect(con) |
|
| 61 |
fs$id=substr(fs$id,7,19) |
|
| 62 |
|
|
| 63 |
### Identify which swaths are available in the datapool |
|
| 64 |
swaths=data.frame(path=list.files(datadir,pattern=paste("hdf$"),recursive=T,full=T),stringsAsFactors=F) #all swaths in data pool
|
|
| 65 |
swaths$id=substr(basename(swaths$path),10,22) |
|
| 66 |
fs$exists=fs$id%in%swaths$id |
|
| 67 |
fs$path=swaths$path[match(fs$id,swaths$id)] |
|
| 68 |
|
|
| 51 |
) |
|
| 52 |
d$tile=tile_bb$tile[i] |
|
| 53 |
print(paste("Finished tile",tile_bb$tile[i]))
|
|
| 54 |
return(d) |
|
| 55 |
})) |
|
| 56 |
con=dbDisconnect(con) |
|
| 57 |
fs$id=substr(fs$id,7,19) |
|
| 58 |
|
|
| 59 |
## Identify which swaths are available in the datapool |
|
| 60 |
swaths=data.frame(path=list.files(datadir,pattern=paste("hdf$"),recursive=T,full=T),stringsAsFactors=F) #all swaths in data pool
|
|
| 61 |
swaths$id=substr(basename(swaths$path),10,22) |
|
| 62 |
fs$exists=fs$id%in%swaths$id |
|
| 63 |
fs$path=swaths$path[match(fs$id,swaths$id)] |
|
| 64 |
|
|
| 65 |
## write tile-swath list to disk |
|
| 66 |
save(fs,swaths,file="swathtile.Rdata") |
|
| 67 |
} |
|
| 68 |
|
|
| 69 |
load("swathtile.Rdata")
|
|
| 70 |
|
|
| 69 | 71 |
if(verbose) print(paste("###############",nrow(fs)," swath IDs recieved from database"))
|
| 70 | 72 |
|
| 71 | 73 |
## get all unique dates |
| 72 | 74 |
fs$dateid=format(as.Date(paste(fs$year,fs$day,sep=""),"%Y%j"),"%Y%m%d") |
| 73 |
alldates=unique(fs$dateid[fs$exists]) |
|
| 75 |
#alldates=unique(fs$dateid[fs$exists])
|
|
| 74 | 76 |
|
| 75 | 77 |
#### Generate submission file |
| 76 | 78 |
startdate="2000-03-01" |
| 77 | 79 |
stopdate="2011-12-31" |
| 78 |
## just 2005 |
|
| 79 |
startdate="2005-01-01"
|
|
| 80 |
stopdate="2005-12-31"
|
|
| 80 |
## just 2005-2010
|
|
| 81 |
startdate="2009-01-01"
|
|
| 82 |
stopdate="2009-12-31"
|
|
| 81 | 83 |
|
| 82 | 84 |
alldates=format(seq(as.Date(startdate),as.Date(stopdate),1),"%Y%m%d") |
| 83 | 85 |
|
| ... | ... | |
| 89 | 91 |
proclist$avail=paste(proclist$tile,proclist$date,sep="_")%in%paste(avail$tile,avail$date,sep="_") |
| 90 | 92 |
|
| 91 | 93 |
## identify which have been completed |
| 92 |
fdone=data.frame(path=system("ssh lou 'find MOD35/daily -name \"*.nc\"' ",intern=T))
|
|
| 93 |
#fdone=data.frame(path=list.files(outdir,pattern="nc$",recursive=T))
|
|
| 94 |
#fdone=data.frame(path=system("ssh lou 'find MOD35/daily -name \"*.nc\"' ",intern=T))
|
|
| 95 |
fdone=data.frame(path=list.files(outdir,pattern="nc$",recursive=T)) |
|
| 94 | 96 |
fdone$date=substr(basename(as.character(fdone$path)),14,21) |
| 95 | 97 |
fdone$tile=substr(basename(as.character(fdone$path)),7,12) |
| 96 | 98 |
proclist$done=paste(proclist$tile,proclist$date,sep="_")%in%substr(basename(as.character(fdone$path)),7,21) |
| ... | ... | |
| 98 | 100 |
### report on what has already been processed |
| 99 | 101 |
print(paste(sum(!proclist$done)," out of ",nrow(proclist)," (",round(100*sum(!proclist$done)/nrow(proclist),2),"%) remain"))
|
| 100 | 102 |
table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done]) |
| 103 |
table(table(tile=proclist$tile[!proclist$done],year=proclist$year[!proclist$done])) |
|
| 104 |
|
|
| 105 |
### explore tile counts |
|
| 106 |
#x=table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done]) |
|
| 107 |
#x=x[order(rownames(x)),] |
|
| 101 | 108 |
|
| 102 | 109 |
script="/u/awilso10/environmental-layers/climate/procedures/MOD35_L2_process.r" |
| 103 | 110 |
|
| ... | ... | |
| 111 | 118 |
### qsub script |
| 112 | 119 |
cat(paste("
|
| 113 | 120 |
#PBS -S /bin/bash |
| 114 |
#PBS -l select=100:ncpus=8:mpiprocs=8 |
|
| 115 |
##PBS -l select=20:ncpus=8:mpiprocs=8
|
|
| 116 |
#PBS -l walltime=5:00:00
|
|
| 117 |
##PBS -l walltime=2:00:00
|
|
| 121 |
##PBS -l select=100:ncpus=8:mpiprocs=8
|
|
| 122 |
#PBS -l select=10:ncpus=8:mpiprocs=8
|
|
| 123 |
##PBS -l walltime=8:00:00
|
|
| 124 |
#PBS -l walltime=2:00:00 |
|
| 118 | 125 |
#PBS -j n |
| 119 | 126 |
#PBS -m be |
| 120 | 127 |
#PBS -N mod35 |
| 121 |
#PBS -q normal |
|
| 122 |
##PBS -q devel
|
|
| 128 |
##PBS -q normal
|
|
| 129 |
#PBS -q devel |
|
| 123 | 130 |
#PBS -V |
| 124 | 131 |
|
| 125 |
CORES=800
|
|
| 132 |
CORES=80 |
|
| 126 | 133 |
#CORES=160 |
| 127 | 134 |
|
| 128 | 135 |
HDIR=/u/armichae/pr/ |
| ... | ... | |
| 147 | 154 |
|
| 148 | 155 |
## Submit it |
| 149 | 156 |
system(paste("qsub mod35_qsub",sep=""))
|
| 157 |
|
|
| 150 | 158 |
system("qstat -u awilso10")
|
| 151 | 159 |
|
| 152 | 160 |
####################################################### |
| 153 | 161 |
### Now submit the script to generate the climatologies |
| 154 | 162 |
|
| 155 | 163 |
tiles |
| 164 |
ctiles=c("h10v08","h11v08","h12v08","h10v07","h11v07","h12v07") # South America
|
|
| 165 |
|
|
| 156 | 166 |
ctiles=tiles#[c(1:3)] #subset to only some tiles (for example if some aren't finished yet)? |
| 157 | 167 |
climatescript="/pleiades/u/awilso10/environmental-layers/climate/procedures/MOD35_Climatology.r" |
| 158 | 168 |
|
| 159 | 169 |
## check which tiles have been processed and are on lou with a filename "MOD35_[tile].nc" |
| 170 |
cdone=data.frame(path="",tile="") #use this if you want to re-run everything |
|
| 160 | 171 |
cdone=data.frame(path=sapply(strsplit(basename( |
| 161 | 172 |
system("ssh lou 'find MOD35/summary -name \"MOD35_h[0-9][0-9]v[0-9][0-9].nc\"' ",intern=T)),split="_"),function(x) x[2]))
|
| 162 | 173 |
cdone$tile=substr(basename(as.character(cdone$path)),1,6) |
| 163 |
print(paste(length(ctiles[!ctiles%in%cdone$tile]),"Tiles still need to be processed: /n ",ctiles[!ctiles%in%cdone$tile]))
|
|
| 174 |
print(paste(length(ctiles[!ctiles%in%cdone$tile]),"Tiles still need to be processed"))
|
|
| 164 | 175 |
|
| 165 | 176 |
## write the table processed by mpiexec |
| 166 | 177 |
write.table(paste("--verbose ",climatescript," --verbose T --tile ",ctiles[!ctiles%in%cdone$tile],sep=""),
|
| ... | ... | |
| 176 | 187 |
### qsub script |
| 177 | 188 |
cat(paste("
|
| 178 | 189 |
#PBS -S /bin/bash |
| 179 |
#PBS -l select=1:ncpus=16:mem=94
|
|
| 180 |
#PBS -l walltime=24:00:00
|
|
| 190 |
#PBS -l select=20:ncpus=8:mem=94
|
|
| 191 |
#PBS -l walltime=3:00:00
|
|
| 181 | 192 |
#PBS -j n |
| 182 | 193 |
#PBS -m be |
| 183 | 194 |
#PBS -N mod35_climate |
| 184 |
#PBS -q ldan |
|
| 195 |
#PBS -q normal |
|
| 196 |
##PBS -q ldan |
|
| 185 | 197 |
#PBS -V |
| 186 | 198 |
",if(delay) paste("#PBS -W depend=afterany:",job,sep=""),"
|
| 187 | 199 |
|
| 188 |
CORES=16 |
|
| 200 |
CORES=160
|
|
| 189 | 201 |
HDIR=/u/armichae/pr/ |
| 190 | 202 |
source $HDIR/etc/environ.sh |
| 191 | 203 |
source /pleiades/u/awilso10/environ.sh |
| ... | ... | |
| 220 | 232 |
|
| 221 | 233 |
system("ssh lou")
|
| 222 | 234 |
#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/ |
| 223 |
rsync -vv `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/
|
|
| 235 |
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/")
|
|
| 224 | 236 |
exit |
| 225 | 237 |
|
| 226 | 238 |
|
Also available in: Unified diff
Move file outputs to /nobackupp1 instead of lou due to quota increase. First full 2009 summary