Revision b6e392da
Added by Adam Wilson about 12 years ago
| climate/procedures/MOD06_L2_process.r | ||
|---|---|---|
| 1 | 1 |
################################################################################### |
| 2 | 2 |
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform |
| 3 | 3 |
|
| 4 |
|
|
| 5 |
# Redirect all warnings to stderr() |
|
| 6 |
#options(warn = -1) |
|
| 7 |
#write("2) write() to stderr", stderr())
|
|
| 8 |
#write("2) write() to stdout", stdout())
|
|
| 9 |
#warning("2) warning()")
|
|
| 10 |
|
|
| 11 |
|
|
| 4 | 12 |
## import commandline arguments |
| 5 | 13 |
library(getopt) |
| 6 | 14 |
## get options |
| 7 | 15 |
opta <- getopt(matrix(c( |
| 8 | 16 |
'date', 'd', 1, 'character', |
| 9 | 17 |
'tile', 't', 1, 'character', |
| 18 |
'verbose','v',1,'logical', |
|
| 10 | 19 |
'help', 'h', 0, 'logical' |
| 11 | 20 |
), ncol=4, byrow=TRUE)) |
| 12 | 21 |
if ( !is.null(opta$help) ) |
| ... | ... | |
| 16 | 25 |
q(status=1); |
| 17 | 26 |
} |
| 18 | 27 |
|
| 19 |
#gridfile = paste("HDF4_EOS:EOS_GRID:\"", opta$data,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl", sep="")
|
|
| 20 |
#cat(paste("Checking : ", gridfile, " ...\n"));
|
|
| 21 |
|
|
| 22 |
|
|
| 23 |
## Load command line arguments (mname) |
|
| 24 |
#args=(commandArgs(TRUE)) ##args is now a list of character vectors |
|
| 25 |
## Then cycle through each element of the list and evaluate the expressions. |
|
| 26 |
#eval(parse(text=args)) |
|
| 27 |
#print(args) |
|
| 28 |
|
|
| 29 |
#system("module list")
|
|
| 30 |
#system("source ~/moduleload")
|
|
| 31 |
#system("module list")
|
|
| 32 |
|
|
| 33 | 28 |
|
| 34 | 29 |
date=opta$date #date="20030301" |
| 35 | 30 |
tile=opta$tile #tile="h11v08" |
| 36 |
outdir="2_daily" #directory for separate daily files
|
|
| 37 |
#outdir2="3_summary" #directory for combined daily files and summarized files
|
|
| 31 |
verbose=opta$verbose #print out extensive information for debugging?
|
|
| 32 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
|
| 38 | 33 |
|
| 39 |
print(paste("Processing tile",tile," for date",date))
|
|
| 34 |
## permanent storage on lou: |
|
| 35 |
ldir=paste("lfe1.nas.nasa.gov:/u/awilso10/mod06/daily/",tile,sep="")
|
|
| 40 | 36 |
|
| 41 |
#system("module list")
|
|
| 42 |
#system("ldd /u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/rgdal/libs/rgdal.so")
|
|
| 37 |
print(paste("Processing tile",tile," for date",date))
|
|
| 43 | 38 |
|
| 44 | 39 |
## load libraries |
| 45 | 40 |
require(reshape) |
| ... | ... | |
| 47 | 42 |
require(raster) |
| 48 | 43 |
library(rgdal) |
| 49 | 44 |
require(spgrass6) |
| 45 |
require(RSQLite) |
|
| 50 | 46 |
|
| 51 | 47 |
|
| 52 | 48 |
## specify some working directories |
| 53 | 49 |
setwd("/nobackupp1/awilso10/mod06")
|
| 54 | 50 |
|
| 55 |
print(paste("tempdir()=",tempdir()))
|
|
| 56 |
print(paste("TMPDIR=",Sys.getenv("TMPDIR")))
|
|
| 57 |
|
|
| 58 | 51 |
## load ancillary data |
| 59 | 52 |
load(file="allfiles.Rdata") |
| 60 | 53 |
|
| 61 | 54 |
## load tile information |
| 62 | 55 |
load(file="modlandTiles.Rdata") |
| 63 |
### use MODIS tile as ROI |
|
| 64 |
#modt=readOGR("modgrid","modis_sinusoidal_grid_world",)
|
|
| 65 |
#modt@data[,colnames(tb)[3:6]]=tb[match(paste(modt$h,modt$v),paste(tb$ih,tb$iv)),3:6] |
|
| 66 |
#write.csv(modt@data,file="modistile.csv") |
|
| 67 | 56 |
|
| 68 | 57 |
## vars to process |
| 69 | 58 |
vars=as.data.frame(matrix(c( |
| ... | ... | |
| 84 | 73 |
############################################################################ |
| 85 | 74 |
### Define functions to process a particular date-tile |
| 86 | 75 |
|
| 76 |
#### get swaths |
|
| 77 |
#getswaths<-function(tile,db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"){
|
|
| 78 |
# db=dbConnect( |
|
| 79 |
# } |
|
| 80 |
|
|
| 87 | 81 |
swath2grid=function(file,vars,upleft,lowright){
|
| 88 | 82 |
## Function to generate hegtool parameter file for multi-band HDF-EOS file |
| 89 | 83 |
print(paste("Starting file",basename(file)))
|
| ... | ... | |
| 118 | 112 |
## write it to a file |
| 119 | 113 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
| 120 | 114 |
## now run the swath2grid tool |
| 121 |
## write the tiff file
|
|
| 122 |
# log=system(paste("/nobackupp4/pvotava/software/heg/bin/swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
|
| 123 |
log=system(paste("/nobackupp1/awilso10/software/heg/bin/swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
|
| 115 |
## write the gridded file and save the log including the pid of the parent process
|
|
| 116 |
# log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T)
|
|
| 117 |
log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
|
| 124 | 118 |
## clean up temporary files in working directory |
| 125 |
# file.remove(paste("filetable.temp_",pid,sep=""))
|
|
| 126 |
print(log) |
|
| 127 |
## confirm file is present |
|
| 128 |
print(paste("Confirming output file (",outfile,") is present and readable by GDAL"))
|
|
| 129 |
system(paste("gdalinfo ",outfile))
|
|
| 119 |
file.remove(list.files(pattern= |
|
| 120 |
paste("filetable.temp_",
|
|
| 121 |
as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process |
|
| 122 |
if(verbose) print(log) |
|
| 130 | 123 |
print(paste("Finished ", file))
|
| 131 | 124 |
} |
| 132 | 125 |
|
| ... | ... | |
| 146 | 139 |
Sys.setenv(GRASS_GUI="txt") |
| 147 | 140 |
|
| 148 | 141 |
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations |
| 149 |
loadcloud<-function(date,fs){
|
|
| 142 |
loadcloud<-function(date,fs,ncfile){
|
|
| 150 | 143 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") |
| 151 |
dir.create(tf) |
|
| 144 |
if(!file.exists(tf))dir.create(tf)
|
|
| 152 | 145 |
|
| 153 | 146 |
print(paste("Set up temporary grass session in",tf))
|
| 154 | 147 |
|
| 155 | 148 |
## set up temporary grass instance for this PID |
| 156 |
# initGRASS(gisBase="/nobackupp1/awilso10/software/grass-6.4.3svn",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
|
| 157 | 149 |
initGRASS(gisBase="/u/armichae/pr/grass-6.4.2/",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
| 158 |
|
|
| 159 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""))
|
|
| 150 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
|
| 160 | 151 |
|
| 161 | 152 |
## Define region by importing one MOD11A1 raster. |
| 162 | 153 |
print("Import one MOD11A1 raster to define grid")
|
| 163 | 154 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
| 164 | 155 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
| 165 |
system("g.region rast=modisgrid save=roi --overwrite")
|
|
| 156 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=T,ignore.stderr=T)
|
|
| 166 | 157 |
|
| 167 | 158 |
## Identify which files to process |
| 168 | 159 |
tfs=fs$file[fs$dateid==date] |
| ... | ... | |
| 177 | 168 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
| 178 | 169 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 179 | 170 |
## extract cloudy and 'confidently clear' pixels |
| 180 |
|
|
| 181 | 171 |
system(paste("r.mapcalc <<EOF
|
| 182 | 172 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
| 183 | 173 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3 |
| 184 | 174 |
EOF",sep="")) |
| 185 |
|
|
| 186 | 175 |
## QA |
| 187 | 176 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| 188 | 177 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| ... | ... | |
| 243 | 232 |
CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),"))*100)
|
| 244 | 233 |
EOF",sep="")) |
| 245 | 234 |
|
| 246 |
#### Write the files to a geotiff |
|
| 247 |
# execGRASS("r.out.gdal",input="CER_daily",output=paste(outdir,"/CER_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
| 248 |
# execGRASS("r.out.gdal",input="COT_daily",output=paste(outdir,"/COT_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
| 249 |
# execGRASS("r.out.gdal",input="CLD_daily",output=paste(outdir,"/CLD_",date,".tif",sep=""),nodata=99,flags=c("quiet"))
|
|
| 250 | 235 |
|
| 251 | 236 |
### Write the files to a netcdf file |
| 252 | 237 |
## create image group to facilitate export as multiband netcdf |
| 253 | 238 |
execGRASS("i.group",group="mod06",input=c("CER_daily","COT_daily","CLD_daily")) ; print("")
|
| 254 | 239 |
|
| 255 |
ncfile=paste(outdir,"/MOD06_",date,".nc",sep="") |
|
| 256 |
file.remove(ncfile) |
|
| 240 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it |
|
| 257 | 241 |
execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),
|
| 258 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
|
| 242 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
|
| 259 | 243 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 260 | 244 |
|
| 261 | 245 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
| 262 | 246 |
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep="")) |
| 263 |
# system(paste("cdo setdate,",as.integer(fs$date[fs$dateid==date]-as.Date("2000-01-01"))," -settaxis,2000-01-01,12:00:00,1day ",ncfile," ",ncfile,sep=""))
|
|
| 264 |
# timdimf=paste(tempdir(),"/time.nc",sep="") |
|
| 265 |
# nc_dtime=ncdim_def("time","days since 2000-01-01 12:00:00",
|
|
| 266 |
# vals=as.integer(fs$date[fs$dateid==date]-as.Date("2000-01-01")),
|
|
| 267 |
# calendar="gregorian",longname="time") |
|
| 268 |
# nc_vtime=ncvar_def("time","
|
|
| 269 |
# nc=nc_create(timdimf,nc_time) |
|
| 270 | 247 |
## create temporary nc file with time information to append to MOD06 data |
| 271 | 248 |
cat(paste("
|
| 272 | 249 |
netcdf time {
|
| ... | ... | |
| 274 | 251 |
time = 1 ; |
| 275 | 252 |
variables: |
| 276 | 253 |
int time(time) ; |
| 277 |
time:units = \"days since 2000-01-01 12:00:00\" ;
|
|
| 254 |
time:units = \"days since 2000-01-01 00:00:00\" ;
|
|
| 278 | 255 |
time:calendar = \"gregorian\"; |
| 279 | 256 |
time:long_name = \"time of observation\"; |
| 280 | 257 |
data: |
| ... | ... | |
| 282 | 259 |
}"),file=paste(tempdir(),"/time.cdl",sep="")) |
| 283 | 260 |
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 284 | 261 |
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep="")) |
| 285 |
|
|
| 286 |
# ncwa -O -h -a record out.nc out.nc |
|
| 287 |
# ncvar_put( nc, "time", vals=as.integer(fs$date[fs$dateid==date]-as.Date("2000-01-01")), start=NA, count=NA, verbose=FALSE )
|
|
| 288 |
# system(paste(ncopath,"ncap2 -O -s 'time[time]=",as.integer(fs$date[fs$dateid==date]-as.Date("2000-01-01")),"'",ncfile," ",ncfile,sep=""))
|
|
| 289 |
# system(paste(ncopath,"ncatted -O calendar,time,c,c,\"standard\" -a long_name,time,c,c,\"time\" -a units,time,c,c,\"days since 2000-01-01 12:00:00\" ",ncfile,sep="")) |
|
| 262 |
## add other attributes |
|
| 290 | 263 |
system(paste(ncopath,"ncrename -v Band1,CER -v Band2,COT -v Band3,CLD ",ncfile,sep="")) |
| 291 | 264 |
system(paste(ncopath,"ncatted -a scale_factor,CER,o,d,0.01 -a units,CER,o,c,\"micron\" -a missing_value,CER,o,d,-32768 -a long_name,CER,o,c,\"Cloud Particle Effective Radius\" ",ncfile,sep="")) |
| 292 | 265 |
system(paste(ncopath,"ncatted -a scale_factor,COT,o,d,0.01 -a units,COT,o,c,\"none\" -a missing_value,COT,o,d,-32768 -a long_name,COT,o,c,\"Cloud Optical Thickness\" ",ncfile,sep="")) |
| ... | ... | |
| 295 | 268 |
|
| 296 | 269 |
### delete the temporary files |
| 297 | 270 |
unlink_.gislock() |
| 298 |
# system("/nobackupp1/awilso10/software/grass-6.4.3svn/etc/clean_temp")
|
|
| 299 |
system(paste("rm -rR ",tf,sep=""))
|
|
| 271 |
system("clean_temp")
|
|
| 272 |
system(paste("rm -frR ",tf,sep=""))
|
|
| 300 | 273 |
} |
| 301 | 274 |
|
| 302 | 275 |
|
| ... | ... | |
| 308 | 281 |
##################################################### |
| 309 | 282 |
## Run the gridding procedure |
| 310 | 283 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
| 284 |
|
|
| 285 |
# if(venezuela) tile_bb$lat_max=11.0780999176 #increase northern boundary for KT |
|
| 286 |
|
|
| 311 | 287 |
lapply(fs$path[fs$dateid==date],swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max)) |
| 312 | 288 |
|
| 289 |
## confirm at least one file for this date is present |
|
| 290 |
outfiles=paste(tempdir(),"/",basename(fs$path[fs$dateid==date]),sep="") |
|
| 291 |
if(!any(file.exists(outfiles))) {
|
|
| 292 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
|
| 293 |
q("no",status=0)
|
|
| 294 |
} |
|
| 295 |
|
|
| 313 | 296 |
##################################################### |
| 314 | 297 |
## Process the gridded files |
| 315 |
|
|
| 316 |
## temporary objects to test function below |
|
| 317 |
# i=1 |
|
| 318 |
#file=paste(outdir,"/",fs$file[1],sep="") |
|
| 319 |
#date=as.Date("2000-05-23")
|
|
| 320 | 298 |
|
| 321 | 299 |
## run themod06 processing for this date |
| 322 |
loadcloud(date,fs=fs) |
|
| 300 |
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
|
| 301 |
loadcloud(date,fs=fs,ncfile=ncfile) |
|
| 302 |
|
|
| 303 |
## Confirm that the file has the correct attributes, otherwise delete it |
|
| 304 |
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
|
| 305 |
npar=as.numeric(system(paste("cdo -s npar ",ncfile),intern=T))
|
|
| 306 |
if(!ntime==1&npar>0) {
|
|
| 307 |
print(paste("File for tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
|
| 308 |
file.remove(ncfile) |
|
| 309 |
} |
|
| 310 |
|
|
| 311 |
## push file to lou and then delete it from pfe |
|
| 312 |
# system(paste("scp ",ncfile," ",ldir,sep="") )
|
|
| 313 |
# file.remove(ncfile) |
|
| 314 |
|
|
| 323 | 315 |
## print out some info |
| 324 | 316 |
print(paste(" ################################################################### Finished ",date,"
|
| 325 | 317 |
################################################################")) |
| ... | ... | |
| 329 | 321 |
##date=notdone[1] |
| 330 | 322 |
mod06(date,tile) |
| 331 | 323 |
|
| 332 |
## run it for all dates |
|
| 324 |
## run it for all dates - Use this if running on a workstation/server (otherwise use qsub)
|
|
| 333 | 325 |
#mclapply(notdone,mod06,tile,mc.cores=ncores) # use ncores/2 because system() commands can add second process for each spawned R |
| 334 | 326 |
#foreach(i=notdone[1:3],.packages=(.packages())) %dopar% mod06(i,tile) |
| 335 | 327 |
#foreach(i=1:20) %dopar% print(i) |
| 336 | 328 |
|
| 329 |
## delete old files |
|
| 330 |
#system("cleartemp")
|
|
| 331 |
|
|
| 337 | 332 |
q("no",status=0)
|
| climate/procedures/Pleiades.R | ||
|---|---|---|
| 1 | 1 |
#### Script to facilitate processing of MOD06 data |
| 2 | 2 |
|
| 3 | 3 |
setwd("/nobackupp1/awilso10/mod06")
|
| 4 |
|
|
| 4 | 5 |
library(rgdal) |
| 5 | 6 |
library(raster) |
| 6 | 7 |
|
| ... | ... | |
| 9 | 10 |
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
|
| 10 | 11 |
save(tb,file="modlandTiles.Rdata") |
| 11 | 12 |
|
| 13 |
## delete temporary log file that can grow to GB |
|
| 14 |
system("rm /nobackupp1/awilso10/software/heg/TOOLKIT_MTD/runtime/LogStatus")
|
|
| 15 |
|
|
| 16 |
tile="h11v08" # Venezuela |
|
| 17 |
tile="h11v07" # Venezuela coast |
|
| 18 |
tile="h09v04" # Oregon |
|
| 19 |
|
|
| 12 | 20 |
outdir="2_daily" #directory for separate daily files |
| 13 |
outdir2="3_summary" #directory for combined daily files and summarized files |
|
| 21 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
|
| 22 |
if(!file.exists(outdir)) dir.create(outdir) |
|
| 23 |
outdir2="summary" #directory for combined daily files and summarized files |
|
| 24 |
if(!file.exists(outdir2)) dir.create(outdir2) |
|
| 14 | 25 |
|
| 15 | 26 |
## load a MOD11A1 file to define grid |
| 16 | 27 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27/",pattern=paste(tile,".*hdf$",sep=""),full=T)[1]
|
| ... | ... | |
| 20 | 31 |
|
| 21 | 32 |
### get list of files to process |
| 22 | 33 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
| 34 |
datadir="/nobackupp1/awilso10/mod06/data" #for data downloaded from |
|
| 23 | 35 |
|
| 24 | 36 |
fs=data.frame(path=list.files(datadir,full=T,recursive=T,pattern="hdf"),stringsAsFactors=F) |
| 25 | 37 |
fs$file=basename(fs$path) |
| ... | ... | |
| 38 | 50 |
|
| 39 | 51 |
#### Generate submission file |
| 40 | 52 |
## identify which have been completed |
| 41 |
done=alldates%in%substr(list.files(outdir),7,14)
|
|
| 42 |
table(done)
|
|
| 43 |
notdone=alldates[!done] #these are the dates that still need to be processed
|
|
| 53 |
#fdone=system(paste("ssh -q lou2 \"ls ",pdir," | grep \"nc$\"\"",sep=""),intern=T)
|
|
| 54 |
fdone=list.files(outdir,pattern="nc$")
|
|
| 55 |
done=alldates%in%substr(fdone,14,21)
|
|
| 44 | 56 |
|
| 45 |
if(exists("fdly")){
|
|
| 46 |
notdone=notdone[notdone%in%fdly$dateid[is.na(fdly$npar)]] |
|
| 57 |
if(exists("fdly")){ #update using table from below
|
|
| 58 |
done[alldates%in%fdly$dateid[is.na(fdly$npar)]]=F |
|
| 59 |
} |
|
| 47 | 60 |
|
| 61 |
table(done) |
|
| 62 |
notdone=alldates[!done] #these are the dates that still need to be processed |
|
| 48 | 63 |
|
| 49 |
tile="h11v08" #can move this to submit script if needed |
|
| 50 | 64 |
script="/u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r" |
| 51 |
write.table(paste("--verbose ",script," --date ",notdone," --tile \"",tile,"\"",sep=""),file="notdone.txt",row.names=F,col.names=F,quote=F)
|
|
| 52 |
#write.table(paste("--verbose ",script," date=",notdone[1:30],sep=""),file="notdone.txt",row.names=F,col.names=F,quote=F)
|
|
| 53 |
#write.table(notdone[1:30],file="notdone.txt",row.names=F,col.names=F,quote=F) |
|
| 65 |
|
|
| 66 |
write.table(paste("--verbose ",script," --date ",notdone," --verbose T --tile \"",tile,"\"",sep=""),file=paste(tile,"_notdone.txt",sep=""),row.names=F,col.names=F,quote=F)
|
|
| 54 | 67 |
|
| 55 | 68 |
save(fs,alldates,gridfile,td,file="allfiles.Rdata") |
| 56 | 69 |
|
| 57 |
## run script |
|
| 58 |
#cat(paste("
|
|
| 59 |
##! /bin/bash |
|
| 60 |
#source ~/moduleload |
|
| 61 |
#source ~/.bashrc |
|
| 62 |
#Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r date=$1 |
|
| 63 |
#Rscript --verbose --vanilla rtest |
|
| 64 |
#",sep=""),file="MOD06_process2") |
|
| 65 |
#system("chmod +x MOD06_process2")
|
|
| 66 |
|
|
| 67 |
#cat(paste("
|
|
| 68 |
#library(rgdal) |
|
| 69 |
#GDALinfo |
|
| 70 |
#",sep=""),file="rtest") |
|
| 71 |
|
|
| 72 |
|
|
| 73 |
## Submission script |
|
| 74 |
|
|
| 75 |
#cat(paste("
|
|
| 76 |
##PBS -S /bin/bash |
|
| 77 |
##PBS -l select=2:ncpus=16:model=san |
|
| 78 |
####PBS -l select=4:ncpus=8:model=neh |
|
| 79 |
###PBS -l select=1:ncpus=12:model=wes |
|
| 80 |
######## old: select=48:ncpus=8:mpiprocs=8:model=neh |
|
| 81 |
##PBS -l walltime=2:00:00 |
|
| 82 |
##PBS -j oe |
|
| 83 |
##PBS -m e |
|
| 84 |
##PBS -V |
|
| 85 |
##PBS -q devel |
|
| 86 |
##PBS -o log/log_^array_index^ |
|
| 87 |
##PBS -o log/log_DataCompile.log |
|
| 88 |
##PBS -M adam.wilson@yale.edu |
|
| 89 |
##PBS -N MOD06 |
|
| 90 |
|
|
| 91 |
## cd to working directory |
|
| 92 |
#cd /nobackupp1/awilso10/mod06 |
|
| 93 |
|
|
| 94 |
## set some memory limits |
|
| 95 |
# ulimit -d 1500000 -m 1500000 -v 1500000 #limit memory usage |
|
| 96 |
# source /usr/local/lib/global.profile |
|
| 97 |
# source /u/awilso10/.bashrc |
|
| 98 |
# source /u/awilso10/moduleload |
|
| 99 |
### export a few important variables |
|
| 100 |
# export NNODES=32 |
|
| 101 |
# export R_LIBS=\"/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/\" |
|
| 102 |
### Run the script! |
|
| 103 |
## current version not parallelizing across nodes! |
|
| 104 |
# TMPDIR=$TMPDIR Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r date=20000403 |
|
| 105 |
|
|
| 106 |
#WORKLIST=notdone.txt |
|
| 107 |
##EXE=\"Rscript\" |
|
| 108 |
#EXE="./MOD06_process2" |
|
| 109 |
#LOG=log/log_DataCompile.log |
|
| 110 |
#MQUEUE=/nobackupp4/pvotava/software/share/mqueue-eg/mqueue/mqueue |
|
| 111 |
|
|
| 112 |
#TMPDIR=$TMPDIR mpiexec -np $NNODES $MQUEUE -l $WORKLIST -p $EXE -v -v -v --random-starts 2-4 --work-analyze #> $LOG |
|
| 113 |
#exit 0 |
|
| 114 |
#",sep=""),file="MOD06_process") |
|
| 115 |
|
|
| 116 |
|
|
| 117 |
### simplified qsub script |
|
| 70 |
### qsub script |
|
| 118 | 71 |
cat(paste("
|
| 119 | 72 |
#PBS -S /bin/bash |
| 120 |
#PBS -l select=48:ncpus=8:mpiprocs=8
|
|
| 73 |
#PBS -l select=40:ncpus=8:mpiprocs=8
|
|
| 121 | 74 |
##PBS -l select=2:ncpus=4:mpiprocs=4 |
| 122 |
#PBS -l walltime=02:00:00
|
|
| 75 |
#PBS -l walltime=2:00:00 |
|
| 123 | 76 |
#PBS -j n |
| 124 | 77 |
#PBS -m be |
| 125 | 78 |
#PBS -N mod06 |
| 126 | 79 |
#PBS -q devel |
| 127 | 80 |
#PBS -V |
| 128 | 81 |
|
| 129 |
CORES=384
|
|
| 82 |
CORES=320
|
|
| 130 | 83 |
HDIR=/u/armichae/pr/ |
| 131 | 84 |
source $HDIR/etc/environ.sh |
| 132 | 85 |
source /u/awilso10/.bashrc |
| 133 | 86 |
IDIR=/nobackupp1/awilso10/mod06/ |
| 134 | 87 |
##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt |
| 135 |
WORKLIST=$IDIR/notdone.txt |
|
| 88 |
WORKLIST=$IDIR/",tile,"_notdone.txt
|
|
| 136 | 89 |
EXE=Rscript |
| 137 | 90 |
LOGSTDOUT=$IDIR/log/log.stdout |
| 138 | 91 |
LOGSTDERR=$IDIR/log/log.stderr |
| 139 | 92 |
mpiexec -np $CORES pxargs -a $WORKLIST -p $EXE -v -v -v --work-analyze 1> $LOGSTDOUT 2> $LOGSTDERR |
| 140 |
"),file="mod06_qsub")
|
|
| 93 |
",sep=""),file=paste(tile,"_mod06_qsub",sep=""))
|
|
| 141 | 94 |
|
| 142 | 95 |
|
| 143 | 96 |
### Check the file |
| 144 |
system("cat mod06_qsub")
|
|
| 97 |
system(paste("cat ",tile,"_mod06_qsub",sep=""))
|
|
| 145 | 98 |
#system("cat ~/environmental-layers/climate/procedures/MOD06_L2_process.r")
|
| 146 | 99 |
|
| 147 | 100 |
## check queue status |
| ... | ... | |
| 149 | 102 |
system("/u/scicon/tools/bin/qtop.pl 492352")
|
| 150 | 103 |
|
| 151 | 104 |
## Submit it (and keep the pid)! |
| 152 |
system("qsub mod06_qsub")
|
|
| 153 |
system("/u/scicon/tools/bin/pdsh_gdb -j 568835 -d tmp -s -u awilso10")
|
|
| 105 |
system(paste("qsub ",tile,"_mod06_qsub",sep=""))
|
|
| 154 | 106 |
|
| 155 | 107 |
## work in interactive mode |
| 156 | 108 |
# system("qsub -I -l walltime=2:00:00 -lselect=2:ncpus=16:model=san -q devel")
|
| ... | ... | |
| 161 | 113 |
system(paste("/u/scicon/tools/bin/qps ",568835))
|
| 162 | 114 |
system(paste("qstat -t -x",pid))
|
| 163 | 115 |
|
| 164 |
system("qstat devel ")
|
|
| 165 |
#system("qstat | grep awilso10")
|
|
| 166 | 116 |
|
| 167 | 117 |
#################################### |
| 168 | 118 |
|
| 119 |
## move to permanent storage on lou? |
|
| 120 |
# system(paste("scp -r ",ncfile," ",pdir,sep=""))
|
|
| 121 |
file.remove( |
|
| 122 |
list.files(pattern="filetable[.]temp|GetAttrtemp|core[.]|MCFWrite[.]temp") |
|
| 123 |
) |
|
| 169 | 124 |
|
| 170 | 125 |
################################################################################ |
| 171 | 126 |
## now generate the climatologies |
| 172 | 127 |
fdly=data.frame( |
| 173 | 128 |
path=list.files(outdir,pattern="nc$",full=T), |
| 174 | 129 |
file=list.files(outdir,pattern="nc$")) |
| 175 |
fdly$dateid=substr(fdly$file,7,14)
|
|
| 176 |
fdly$date=as.Date(substr(fdly$file,7,14),"%Y%m%d")
|
|
| 130 |
fdly$dateid=substr(fdly$file,14,21)
|
|
| 131 |
fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")
|
|
| 177 | 132 |
fdly$month=format(fdly$date,"%m") |
| 178 | 133 |
fdly$year=format(fdly$date,"%Y") |
| 134 |
nrow(fdly) |
|
| 179 | 135 |
|
| 180 | 136 |
## check validity (via npar and ntime) of nc files |
| 181 | 137 |
for(i in 1:nrow(fdly)){
|
| 182 |
# fdly$ntime[i]<-as.numeric(system(paste("cdo sinfo ",fdly$path[i]),intern=T))
|
|
| 138 |
fdly$ntime[i]<-as.numeric(system(paste("cdo -s ntime ",fdly$path[i]),intern=T))
|
|
| 183 | 139 |
fdly$npar[i]<-as.numeric(system(paste("cdo -s npar ",fdly$path[i]),intern=T))
|
| 184 |
print(i) |
|
| 140 |
fdly$fyear[i]<-as.numeric(system(paste("cdo -s showyear ",fdly$path[i]),intern=T))
|
|
| 141 |
fdly$fmonth[i]<-as.numeric(system(paste("cdo -s showmon ",fdly$path[i]),intern=T))
|
|
| 142 |
fdly$fvar[i]<-system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
|
|
| 143 |
print(paste(i," out of ",nrow(fdly)," for year ", fdly$fyear[i])) |
|
| 144 |
} |
|
| 145 |
|
|
| 146 |
### table of problematic files |
|
| 147 |
table(is.na(fdly$npar)) |
|
| 148 |
table(fdly$fmonth) |
|
| 149 |
table(fdly$fvar) |
|
| 150 |
|
|
| 151 |
fdly[is.na(fdly$npar),] |
|
| 152 |
|
|
| 153 |
## delete files that fail check? |
|
| 154 |
delete=T |
|
| 155 |
if(delete) {
|
|
| 156 |
drop=is.na(fdly$npar)|fdly$fvar!=" CER COT CLD" |
|
| 157 |
file.remove(as.character(fdly$path[drop])) |
|
| 158 |
fdly=fdly[!drop,] |
|
| 185 | 159 |
} |
| 186 | 160 |
|
| 187 | 161 |
## Combine all days within years into a single file (can't mergetime all days at once because this opens too many files) |
| 162 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
|
| 163 |
|
|
| 164 |
library(multicore) |
|
| 165 |
|
|
| 166 |
## create temporary directory to put intermediate files (will be deleted when R quits) |
|
| 188 | 167 |
tsdir=paste(tempdir(),"/summary",sep="") |
| 189 | 168 |
dir.create(tsdir) |
| 190 |
lapply(unique(fdly$year),function(y){
|
|
| 191 |
system(paste("cdo -O mergetime ",paste(fdly$path[!is.na(fdly$npar)&fdly$year==y],collapse=" ")," ",tsdir,"/MOD09_",tile,"_",y,"_daily.nc",sep=""))
|
|
| 192 |
print(paste("Finished merging daily files for year",y))
|
|
| 193 |
}) |
|
| 194 |
## Combine the year-by-year files into a single daily file |
|
| 195 |
system(paste("cdo -O mergetime ",paste(list.files(tsdir,full=T,pattern="daily[.]nc$"),collapse=" ")," ",outdir2,"/MOD09_",tile,"_daily.nc",sep=""))
|
|
| 196 | 169 |
|
| 197 |
system(paste("cdo -O monmean ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_monmean.nc",sep=""))
|
|
| 198 |
system(paste("cdo -O ymonmean ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_ymonmean.nc",sep=""))
|
|
| 199 |
system(paste("cdo -O ymonstd ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_ymonstd.nc",sep=""))
|
|
| 170 |
## create a file of daily scenes for each year |
|
| 171 |
#mclapply(unique(fdly$year),function(y){
|
|
| 172 |
# system(paste("cdo -O mergetime ",paste(fdly$path[!is.na(fdly$npar)&fdly$year==y],collapse=" ")," ",tsdir,"/MOD06_",tile,"_",y,"_daily.nc",sep=""))
|
|
| 173 |
# print(paste("Finished merging daily files for year",y))
|
|
| 174 |
#},mc.cores=5) |
|
| 175 |
|
|
| 176 |
|
|
| 177 |
## merge all daily files to create a single file with all dates |
|
| 178 |
system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD06_",tile,"_daily.nc",sep="")) |
|
| 179 |
|
|
| 180 |
############################## |
|
| 181 |
## Combine the year-by-year files into a single daily file in the summary directory (for archiving) |
|
| 182 |
#system(paste("cdo -O mergetime ",paste(list.files(tsdir,full=T,pattern="daily[.]nc$"),collapse=" ")," ",outdir2,"/MOD06_",tile,"_daily.nc",sep=""))
|
|
| 183 |
system(paste(ncopath,"ncatted ", |
|
| 184 |
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ", |
|
| 185 |
" -a title,global,o,c,\"MODIS Cloud Product (MOD06) Daily Timeseries\" ", |
|
| 186 |
" -a institution,global,o,c,\"Yale University\" ", |
|
| 187 |
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ", |
|
| 188 |
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ", |
|
| 189 |
outdir2,"/MOD06_",tile,"_daily.nc",sep="")) |
|
| 190 |
|
|
| 191 |
#system(paste("cdo -O monmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_monmean.nc",sep=""))
|
|
| 192 |
|
|
| 193 |
############################# |
|
| 194 |
## Generate the Climatologies |
|
| 195 |
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15) |
|
| 196 |
|
|
| 197 |
## Monthly means |
|
| 198 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""))
|
|
| 199 |
|
|
| 200 |
## Monthly standard deviation |
|
| 201 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
|
|
| 202 |
system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep="")) |
|
| 203 |
system(paste(ncopath,"ncatted ", |
|
| 204 |
" -a long_name,CER_sd,o,c,\"Cloud Particle Effective Radius (standard deviation of daily observations)\" ", |
|
| 205 |
" -a long_name,CLD_sd,o,c,\"Cloud Mask (standard deviation of daily observations)\" ", |
|
| 206 |
" -a long_name,COT_sd,o,c,\"Cloud Optical Thickness (standard deviation of daily observations)\" ", |
|
| 207 |
tsdir,"/MOD06_",tile,"_ymonstd.nc",sep="")) |
|
| 208 |
|
|
| 209 |
## cer > 20 |
|
| 210 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -gtc,20 -selvar,CER ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
|
|
| 211 |
system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep="")) |
|
| 212 |
system(paste(ncopath,"ncatted ", |
|
| 213 |
" -a long_name,CER20,o,c,\"Proportion of Days with Cloud Particle Effective Radius > 20um\" ", |
|
| 214 |
" -a units,CER20,o,c,\"Proportion\" ", |
|
| 215 |
tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep="")) |
|
| 216 |
|
|
| 217 |
## number of observations |
|
| 218 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,9999 -setmisstoc,9999 -selvar,CER,CLD ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
|
|
| 219 |
system(paste(ncopath,"ncrename -v CER,CER_pmiss -v CLD,CLD_pmiss ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep="")) |
|
| 220 |
system(paste(ncopath,"ncatted ", |
|
| 221 |
" -a long_name,CER_pmiss,o,c,\"Proportion of Days with missing data for CER\" ", |
|
| 222 |
" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ", |
|
| 223 |
" -a scale_factor,CER_pmiss,o,d,0.01 ", |
|
| 224 |
" -a units,CER_pmiss,o,c,\"Proportion\" ", |
|
| 225 |
" -a scale_factor,CLD_pmiss,o,d,0.01 ", |
|
| 226 |
" -a units,CLD_pmiss,o,c,\"Proportion\" ", |
|
| 227 |
tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep="")) |
|
| 228 |
|
|
| 229 |
## append variables to a single file |
|
| 230 |
system(paste(ncopath,"ncks -O ",tsdir,"/MOD06_",tile,"_ymonmean.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 231 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonstd.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 232 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncer20.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 233 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 234 |
|
|
| 235 |
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes |
|
| 236 |
system(paste(ncopath,"ncks -A -v sinusoidal ",list.files(tsdir,full=T,pattern="daily[.]nc$")[1]," ",tsdir,"/MOD06_",tile,"_ymon.nc",sep="")) |
|
| 237 |
|
|
| 238 |
## invert latitude so it plays nicely with gdal |
|
| 239 |
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD06_",tile,"_ymon.nc ",outdir2,"/MOD06_",tile,".nc",sep="")) |
|
| 240 |
|
|
| 241 |
## update attributes |
|
| 242 |
system(paste(ncopath,"ncatted ", |
|
| 243 |
#" -a standard_parallel,sinusoidal,o,c,\"0\" ", |
|
| 244 |
#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ", |
|
| 245 |
#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ", |
|
| 246 |
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ", |
|
| 247 |
" -a title,global,o,c,\"MODIS Cloud Product (MOD06) Climatology\" ", |
|
| 248 |
" -a institution,global,o,c,\"Yale University\" ", |
|
| 249 |
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ", |
|
| 250 |
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ", |
|
| 251 |
outdir2,"/MOD06_",tile,".nc",sep="")) |
|
| 252 |
|
|
| 200 | 253 |
|
| 201 | 254 |
print("Finished! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
|
| 202 | 255 |
## quit R |
| ... | ... | |
| 206 | 259 |
################################################################# |
| 207 | 260 |
|
| 208 | 261 |
### copy the files back to Yale |
| 209 |
list.files("2_daily")
|
|
| 210 |
system(paste("scp ",outdir2,"/*_ymonmean.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela/summary",sep=""))
|
|
| 262 |
system(paste("scp ",outdir2,"/MOD06_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/summary",sep=""))
|
|
| 211 | 263 |
|
| 212 | 264 |
system("scp /tmp/Rtmp6I6tFn/MOD06_L2.A2000061.1615.051.2010273184629.hdf adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela")
|
| 213 | 265 |
system("scp 2_daily/MOD06_20000410.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela")
|
Also available in: Unified diff
Identified problem in handling of missing values of MOD06 product. In some areas limiting quality control to the highest quality results in ~90% NAs.