Revision be64daa8
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_Climatology.r | ||
|---|---|---|
| 33 | 33 |
|
| 34 | 34 |
################################################################################ |
| 35 | 35 |
## Get list of all daily files |
| 36 |
if(verbose) print("Checking daily output in preparation for generating climatology")
|
|
| 36 |
if(verbose) print(paste("Checking daily output in preparation for generating climatology:",tile))
|
|
| 37 | 37 |
|
| 38 | 38 |
fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F) |
| 39 | 39 |
fdly$file=basename(fdly$path) |
| ... | ... | |
| 44 | 44 |
nrow(fdly) |
| 45 | 45 |
|
| 46 | 46 |
## print some summaries |
| 47 |
if(verbose) print("Summary of available daily files")
|
|
| 47 |
if(verbose) print(paste("Summary of available daily files:",tile))
|
|
| 48 | 48 |
print(table(fdly$year)) |
| 49 | 49 |
print(table(fdly$month)) |
| 50 | 50 |
#print(table(fdly$fvar)) |
| 51 | 51 |
|
| 52 | 52 |
################################################################################# |
| 53 | 53 |
## Combine the year-by-year files into a single daily file in the summary directory (for archiving) |
| 54 |
if(verbose) print("Merging daily files into single file output")
|
|
| 54 |
if(verbose) print(paste("Merging daily files into single file output:",tile))
|
|
| 55 | 55 |
|
| 56 | 56 |
## create temporary directory to put intermediate files (will be deleted when R quits) |
| 57 | 57 |
tsdir=paste(tempdir(),"/summary",sep="") |
| ... | ... | |
| 87 | 87 |
|
| 88 | 88 |
############################# |
| 89 | 89 |
## Generate the Climatologies |
| 90 |
if(verbose) print("Generate monthly climatologies")
|
|
| 90 |
if(verbose) print(paste("Generate monthly climatologies: ",tile))
|
|
| 91 | 91 |
|
| 92 | 92 |
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15) |
| 93 | 93 |
|
| 94 | 94 |
## Monthly means |
| 95 |
if(verbose) print("Calculating the monthly means")
|
|
| 95 |
if(verbose) print(paste("Calculating the monthly means:",tile))
|
|
| 96 | 96 |
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)
|
| 97 | 97 |
system(paste(ncopath,"ncrename -v PClear,PCloud ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep="")) |
| 98 | 98 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 112 | 112 |
|
| 113 | 113 |
|
| 114 | 114 |
## Monthly standard deviation |
| 115 |
if(verbose) print("Calculating the monthly SD")
|
|
| 116 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -monmean ",
|
|
| 115 |
if(verbose) print(paste("Calculating the monthly SD:",tile))
|
|
| 116 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -mulc,-1 -subc,100 -monmean ",
|
|
| 117 | 117 |
outdir2,"/MOD35_",tile,"_daily.nc ", |
| 118 | 118 |
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep="")) |
| 119 | 119 |
system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep="")) |
| ... | ... | |
| 122 | 122 |
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep="")) |
| 123 | 123 |
|
| 124 | 124 |
## frequency of cloud days p(clear<90%) |
| 125 |
if(verbose) print("Calculating the proportion of cloudy and probably cloudy days")
|
|
| 125 |
if(verbose) print(paste("Calculating the proportion of cloudy and probably cloudy days:",tile))
|
|
| 126 | 126 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -lec,90 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 127 | 127 |
system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep="")) |
| 128 | 128 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 131 | 131 |
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep="")) |
| 132 | 132 |
|
| 133 | 133 |
## number of observations |
| 134 |
if(verbose) print("Calculating the number of missing variables")
|
|
| 134 |
if(verbose) print(paste("Calculating the number of missing variables:",tile))
|
|
| 135 | 135 |
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -eqc,9999 -setmisstoc,9999 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 136 | 136 |
system(paste(ncopath,"ncrename -v PClear,Pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep="")) |
| 137 | 137 |
system(paste(ncopath,"ncatted ", |
| ... | ... | |
| 143 | 143 |
## clean up variables? |
| 144 | 144 |
|
| 145 | 145 |
## append variables to a single file |
| 146 |
if(verbose) print("Append all monthly climatologies into a single file")
|
|
| 146 |
if(verbose) print(paste("Append all monthly climatologies into a single file:",tile))
|
|
| 147 | 147 |
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep="")) |
| 148 | 148 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep="")) |
| 149 | 149 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep="")) |
| 150 | 150 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep="")) |
| 151 | 151 |
|
| 152 | 152 |
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes |
| 153 |
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")
|
|
| 153 |
if(verbose) print(paste("Clean up file (update attributes, flip latitudes, add grid description:",tile))
|
|
| 154 | 154 |
|
| 155 | 155 |
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep="")) |
| 156 | 156 |
#system(paste(ncopath,"ncks -A -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/MOD35_",tile,"_ymon.nc",sep="")) |
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 93 | 93 |
## load tile information and get bounding box |
| 94 | 94 |
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata") |
| 95 | 95 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
| 96 |
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner |
|
| 97 |
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner |
|
| 96 |
|
|
| 97 |
## get bounds of swath to keep and feed into grass when generating tile |
|
| 98 |
## expand a little (0.5 deg) to ensure that there is no clipping of pixels on the edges |
|
| 99 |
## tile will later be aligned with MODLAND tile so the extra will eventually be trimmed |
|
| 100 |
upleft=paste(min(90,tile_bb$lat_max+.5),max(-180,tile_bb$lon_min-.5)) #northwest corner |
|
| 101 |
lowright=paste(max(-90,tile_bb$lat_min-0.5),min(180,tile_bb$lon_max+0.5)) #southeast corner |
|
| 98 | 102 |
|
| 99 | 103 |
|
| 100 | 104 |
## vars to process |
| climate/procedures/Pleiades_MOD35.R | ||
|---|---|---|
| 75 | 75 |
#### Generate submission file |
| 76 | 76 |
startdate="2000-03-01" |
| 77 | 77 |
stopdate="2011-12-31" |
| 78 |
## just 2009
|
|
| 79 |
startdate="2009-01-01"
|
|
| 80 |
stopdate="2009-12-31"
|
|
| 78 |
## just 2005
|
|
| 79 |
startdate="2005-01-01"
|
|
| 80 |
stopdate="2005-12-31"
|
|
| 81 | 81 |
|
| 82 | 82 |
alldates=format(seq(as.Date(startdate),as.Date(stopdate),1),"%Y%m%d") |
| 83 | 83 |
|
| ... | ... | |
| 102 | 102 |
script="/u/awilso10/environmental-layers/climate/procedures/MOD35_L2_process.r" |
| 103 | 103 |
|
| 104 | 104 |
## write the table processed by mpiexec |
| 105 |
tp=(!proclist$done)&proclist$avail #date-tiles to process
|
|
| 105 |
tp=((!proclist$done)&proclist$avail) #date-tiles to process
|
|
| 106 | 106 |
table(Available=proclist$avail,Completed=proclist$done) |
| 107 | 107 |
|
| 108 | 108 |
write.table(paste("--verbose ",script," --date ",proclist$date[tp]," --verbose T --tile ",proclist$tile[tp],sep=""),
|
| ... | ... | |
| 111 | 111 |
### qsub script |
| 112 | 112 |
cat(paste("
|
| 113 | 113 |
#PBS -S /bin/bash |
| 114 |
#PBS -l select=50:ncpus=8:mpiprocs=8 |
|
| 115 |
##PBS -l select=2:ncpus=8:mpiprocs=8 |
|
| 116 |
##PBS -l select=2:ncpus=4:mpiprocs=4 |
|
| 114 |
#PBS -l select=100:ncpus=8:mpiprocs=8 |
|
| 115 |
##PBS -l select=20:ncpus=8:mpiprocs=8 |
|
| 117 | 116 |
#PBS -l walltime=5:00:00 |
| 117 |
##PBS -l walltime=2:00:00 |
|
| 118 | 118 |
#PBS -j n |
| 119 | 119 |
#PBS -m be |
| 120 | 120 |
#PBS -N mod35 |
| 121 | 121 |
#PBS -q normal |
| 122 |
##PBS -q devel |
|
| 122 | 123 |
#PBS -V |
| 123 | 124 |
|
| 124 |
CORES=400 |
|
| 125 |
CORES=800 |
|
| 126 |
#CORES=160 |
|
| 127 |
|
|
| 125 | 128 |
HDIR=/u/armichae/pr/ |
| 126 | 129 |
# source $HDIR/etc/environ.sh |
| 127 | 130 |
source /u/awilso10/environ.sh |
| ... | ... | |
| 150 | 153 |
### Now submit the script to generate the climatologies |
| 151 | 154 |
|
| 152 | 155 |
tiles |
| 153 |
ctiles=tiles[c(1:3)] #subset to only some tiles (for example if some aren't finished yet)? |
|
| 156 |
ctiles=tiles#[c(1:3)] #subset to only some tiles (for example if some aren't finished yet)?
|
|
| 154 | 157 |
climatescript="/pleiades/u/awilso10/environmental-layers/climate/procedures/MOD35_Climatology.r" |
| 155 | 158 |
|
| 159 |
## check which tiles have been processed and are on lou with a filename "MOD35_[tile].nc" |
|
| 160 |
cdone=data.frame(path=sapply(strsplit(basename( |
|
| 161 |
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 |
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])) |
|
| 164 |
|
|
| 156 | 165 |
## write the table processed by mpiexec |
| 157 |
write.table(paste("--verbose ",climatescript," --verbose T --tile ",ctiles,sep=""),
|
|
| 166 |
write.table(paste("--verbose ",climatescript," --verbose T --tile ",ctiles[!ctiles%in%cdone$tile],sep=""),
|
|
| 158 | 167 |
file=paste("notdone_climate.txt",sep=""),row.names=F,col.names=F,quote=F)
|
| 159 | 168 |
|
| 160 | 169 |
## delay start until previous jobs have finished? |
| ... | ... | |
| 177 | 186 |
",if(delay) paste("#PBS -W depend=afterany:",job,sep=""),"
|
| 178 | 187 |
|
| 179 | 188 |
CORES=16 |
| 180 |
HDIR=/pleiades/u/armichae/pr/
|
|
| 189 |
HDIR=/u/armichae/pr/ |
|
| 181 | 190 |
source $HDIR/etc/environ.sh |
| 182 | 191 |
source /pleiades/u/awilso10/environ.sh |
| 183 | 192 |
source /pleiades/u/awilso10/.bashrc |
| 184 |
source /pleiades/u/awilso10/moduleload |
|
| 185 | 193 |
IDIR=/nobackupp1/awilso10/mod35/ |
| 186 | 194 |
##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt |
| 187 | 195 |
WORKLIST=$IDIR/notdone_climate.txt |
| ... | ... | |
| 209 | 217 |
|
| 210 | 218 |
################################################################# |
| 211 | 219 |
### copy the files back to Yale |
| 212 |
summarydir="summary" |
|
| 213 |
|
|
| 214 |
sumfiles=list.files("summary",pattern="^MOD06_.*[0-9][.]nc",full=T)
|
|
| 215 |
|
|
| 216 |
system(paste("scp ",paste(sumfiles,collapse=" ")," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/summary",sep=""))
|
|
| 217 |
|
|
| 218 |
#system(paste("scp ",tsdir,"/MOD06_",tile,"*.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/summary",sep=""))
|
|
| 219 |
#system(paste("scp ",paste(fs$path[40421:40422],collapse=" ")," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/swaths",sep=""))
|
|
| 220 |
|
|
| 221 |
|
|
| 222 | 220 |
|
| 221 |
system("ssh lou")
|
|
| 222 |
#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/ |
|
| 224 |
exit |
|
| 223 | 225 |
|
| 224 | 226 |
|
Also available in: Unified diff
Updated swath-grid section of MOD35 processing to add buffer and (hopefully) reduce edge artifacts in tiles