Revision 3eb8b23f
Added by Adam Wilson almost 12 years ago
| climate/procedures/MOD06_L2_process.r | ||
|---|---|---|
| 26 | 26 |
} |
| 27 | 27 |
|
| 28 | 28 |
|
| 29 |
date=opta$date #date="20030301" |
|
| 30 |
tile=opta$tile #tile="h11v08" |
|
| 29 |
## default date and tile to play with |
|
| 30 |
date="20030301" |
|
| 31 |
tile="h11v08" |
|
| 32 |
platform="pleiades" |
|
| 33 |
verbose=T |
|
| 34 |
|
|
| 35 |
## now update using options if given |
|
| 36 |
date=opta$date |
|
| 37 |
tile=opta$tile |
|
| 31 | 38 |
verbose=opta$verbose #print out extensive information for debugging? |
| 32 | 39 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 33 | 40 |
|
| 34 |
## location of MOD06 files |
|
| 35 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 41 |
|
|
| 42 |
if(platform="pleiades"){
|
|
| 43 |
## location of MOD06 files |
|
| 44 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 45 |
## path to some executables |
|
| 46 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
|
| 47 |
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif" |
|
| 48 |
## path to swath database |
|
| 49 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 50 |
## specify working directory |
|
| 51 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 52 |
gisBase="/u/armichae/pr/grass-6.4.2/" |
|
| 53 |
## path to MOD11A1 file for this tile to align grid/extent |
|
| 54 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
|
| 55 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
|
| 56 |
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 " |
|
| 57 |
} |
|
| 58 |
|
|
| 59 |
if(platform="litoria"){ #if running on local server, use different paths
|
|
| 60 |
## location of MOD06 files |
|
| 61 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 62 |
## path to some executables |
|
| 63 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
|
| 64 |
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif" |
|
| 65 |
## path to swath database |
|
| 66 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 67 |
## specify working directory |
|
| 68 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 69 |
gisBase="/usr/lib/grass64" |
|
| 70 |
## get grid file |
|
| 71 |
td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
|
| 72 |
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 " |
|
| 73 |
} |
|
| 74 |
|
|
| 36 | 75 |
|
| 37 | 76 |
### print some status messages |
| 38 |
print(paste("Processing tile",tile," for date",date))
|
|
| 77 |
if(verbose) print(paste("Processing tile",tile," for date",date))
|
|
| 39 | 78 |
|
| 40 | 79 |
## load libraries |
| 41 | 80 |
require(reshape) |
| ... | ... | |
| 45 | 84 |
require(spgrass6) |
| 46 | 85 |
require(RSQLite) |
| 47 | 86 |
|
| 48 |
## path to swath database |
|
| 49 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 50 |
## path to NCO |
|
| 51 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
|
| 52 |
|
|
| 53 |
## specify working directory |
|
| 54 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 55 | 87 |
|
| 56 |
## load tile information |
|
| 88 |
## load tile information and get bounding box
|
|
| 57 | 89 |
load(file="modlandTiles.Rdata") |
| 90 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
|
| 91 |
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner |
|
| 92 |
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner |
|
| 58 | 93 |
|
| 59 |
## path to MOD11A1 file for this tile to align grid/extent |
|
| 60 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
|
| 61 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
|
| 62 |
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 " |
|
| 63 | 94 |
|
| 64 | 95 |
## vars to process |
| 65 | 96 |
vars=as.data.frame(matrix(c( |
| ... | ... | |
| 78 | 109 |
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end. |
| 79 | 110 |
finalvars=c("CER","COT","CLD")
|
| 80 | 111 |
|
| 112 |
|
|
| 113 |
##################################################### |
|
| 114 |
##find swaths in region from sqlite database for the specified date/tile |
|
| 115 |
if(verbose) print("Accessing swath ID's from database")
|
|
| 116 |
con=dbConnect("SQLite", dbname = db)
|
|
| 117 |
fs=dbGetQuery(con,paste("SELECT * from swath_geo
|
|
| 118 |
WHERE east>=",tile_bb$lon_min," AND |
|
| 119 |
west<=",tile_bb$lon_max," AND |
|
| 120 |
north>=",tile_bb$lat_min," AND |
|
| 121 |
south<=",tile_bb$lat_max," AND |
|
| 122 |
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND |
|
| 123 |
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j")) |
|
| 124 |
)) |
|
| 125 |
con=dbDisconnect(con) |
|
| 126 |
fs$id=substr(fs$id,7,19) |
|
| 127 |
## find the swaths on disk (using datadir) |
|
| 128 |
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T) |
|
| 129 |
|
|
| 130 |
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk")) |
|
| 131 |
|
|
| 132 |
|
|
| 81 | 133 |
############################################################################ |
| 82 | 134 |
############################################################################ |
| 83 |
### Define functions to process a particular date-tile |
|
| 84 |
|
|
| 85 |
swath2grid=function(file,vars,upleft,lowright){
|
|
| 135 |
### Use the HEG tool to grid all available swath data for this date-tile |
|
| 136 |
for(file in swaths){
|
|
| 86 | 137 |
## Function to generate hegtool parameter file for multi-band HDF-EOS file |
| 87 | 138 |
print(paste("Starting file",basename(file)))
|
| 88 | 139 |
outfile=paste(tempdir(),"/",basename(file),sep="") |
| ... | ... | |
| 116 | 167 |
## write it to a file |
| 117 | 168 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
| 118 | 169 |
## now run the swath2grid tool |
| 119 |
## write the gridded file and save the log including the pid of the parent process |
|
| 120 |
# log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T)
|
|
| 121 |
log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
|
| 170 |
## write the gridded file |
|
| 171 |
log=system(paste("(",swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
|
| 122 | 172 |
## clean up temporary files in working directory |
| 123 |
file.remove(list.files(pattern= |
|
| 124 |
paste("filetable.temp_",
|
|
| 125 |
as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process |
|
| 173 |
# file.remove(list.files(pattern=
|
|
| 174 |
# paste("filetable.temp_",
|
|
| 175 |
# as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process
|
|
| 126 | 176 |
if(verbose) print(log) |
| 127 |
print(paste("Finished ", file))
|
|
| 177 |
print(paste("Finished gridding ", file))
|
|
| 178 |
} #end looping over swaths |
|
| 179 |
|
|
| 180 |
######################## |
|
| 181 |
## confirm at least one file for this date is present. If not, quit. |
|
| 182 |
outfiles=paste(tempdir(),"/",basename(swaths),sep="") |
|
| 183 |
if(!any(file.exists(outfiles))) {
|
|
| 184 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
|
| 185 |
q("no",status=0)
|
|
| 128 | 186 |
} |
| 129 | 187 |
|
| 130 |
|
|
| 131 |
############################################################## |
|
| 132 |
### Import to GRASS for processing |
|
| 188 |
##################################################### |
|
| 189 |
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths |
|
| 190 |
|
|
| 191 |
## Identify output file |
|
| 192 |
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
|
| 133 | 193 |
|
| 134 | 194 |
## function to convert binary to decimal to assist in identifying correct values |
| 135 | 195 |
## this is helpful when defining QA handling below, but isn't used in processing |
| ... | ... | |
| 142 | 202 |
Sys.setenv(DEBUG=1) |
| 143 | 203 |
Sys.setenv(GRASS_GUI="txt") |
| 144 | 204 |
|
| 145 |
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations |
|
| 146 |
loadcloud<-function(date,swaths,ncfile){
|
|
| 147 |
## make temporary working directory |
|
| 205 |
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations |
|
| 206 |
## make temporary working directory |
|
| 148 | 207 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar |
| 149 | 208 |
if(!file.exists(tf)) dir.create(tf) |
| 150 | 209 |
## create output directory if needed |
| ... | ... | |
| 152 | 211 |
|
| 153 | 212 |
## set up temporary grass instance for this PID |
| 154 | 213 |
if(verbose) print(paste("Set up temporary grass session in",tf))
|
| 155 |
initGRASS(gisBase="/u/armichae/pr/grass-6.4.2/",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 214 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 156 | 215 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 157 | 216 |
|
| 158 | 217 |
## Define region by importing one MOD11A1 raster. |
| 159 | 218 |
print("Import one MOD11A1 raster to define grid")
|
| 160 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 219 |
if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 161 | 220 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
| 162 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=T,ignore.stderr=T)
|
|
| 221 |
if(platform=="litoria") execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 222 |
|
|
| 223 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=T,ignore.stderr=T)
|
|
| 163 | 224 |
|
| 164 |
## Identify which files to process |
|
| 165 |
tfs=basename(swaths) |
|
| 166 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
|
| 167 |
tfs=tfs[tfs%in%list.files(tempdir())] |
|
| 168 |
nfs=length(tfs) |
|
| 225 |
## Identify which files to process |
|
| 226 |
tfs=basename(swaths) |
|
| 227 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
|
| 228 |
tfs=tfs[tfs%in%list.files(tempdir())] |
|
| 229 |
nfs=length(tfs) |
|
| 230 |
if(verbose) print(nfs,"swaths available for processing") |
|
| 169 | 231 |
|
| 170 |
## loop through scenes and process QA flags
|
|
| 232 |
## loop through scenes and process QA flags |
|
| 171 | 233 |
for(i in 1:nfs){
|
| 172 | 234 |
file=paste(tempdir(),"/",tfs[i],sep="") |
| 173 | 235 |
## Cloud Mask |
| ... | ... | |
| 178 | 240 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
| 179 | 241 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) > 2 |
| 180 | 242 |
EOF",sep="")) |
| 181 |
## QA |
|
| 243 |
|
|
| 244 |
## QA |
|
| 182 | 245 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| 183 | 246 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 184 | 247 |
## QA_CER |
| ... | ... | |
| 199 | 262 |
flags=c("overwrite","o")) ; print("")
|
| 200 | 263 |
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
| 201 | 264 |
## keep only positive COT values where quality is 'useful' and '>= good' & scale to real units |
| 202 |
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 265 |
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,",null())\"",sep=""))
|
|
| 203 | 266 |
## set COT to 0 in clear-sky pixels |
| 204 | 267 |
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
| 205 | 268 |
|
| ... | ... | |
| 210 | 273 |
flags=c("overwrite","o")) ; print("")
|
| 211 | 274 |
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
| 212 | 275 |
## keep only positive CER values where quality is 'useful' and '>= good' & scale to real units |
| 213 |
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 276 |
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,",null())\"",sep=""))
|
|
| 214 | 277 |
## set CER to 0 in clear-sky pixels |
| 215 | 278 |
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
| 216 | 279 |
|
| ... | ... | |
| 220 | 283 |
# flags=c("overwrite","o")) ; print("")
|
| 221 | 284 |
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
| 222 | 285 |
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units |
| 223 |
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 286 |
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,",null())\"",sep=""))
|
|
| 224 | 287 |
## set CER to 0 in clear-sky pixels |
| 225 | 288 |
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
| 226 | 289 |
|
| ... | ... | |
| 231 | 294 |
system(paste("r.mapcalc <<EOF
|
| 232 | 295 |
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
| 233 | 296 |
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
| 234 |
COT_daily=int((COT_numer/COT_denom)*100)
|
|
| 297 |
COT_daily=int(COT_numer/COT_denom)
|
|
| 235 | 298 |
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
| 236 | 299 |
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
| 237 |
CER_daily=int(100*(CER_numer/CER_denom))
|
|
| 300 |
CER_daily=int(CER_numer/CER_denom)
|
|
| 238 | 301 |
CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),"))*100)
|
| 239 | 302 |
EOF",sep="")) |
| 240 | 303 |
|
| ... | ... | |
| 248 | 311 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 249 | 312 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 250 | 313 |
|
| 251 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
|
| 252 | 314 |
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep="")) |
| 253 | 315 |
## create temporary nc file with time information to append to MOD06 data |
| 254 | 316 |
cat(paste("
|
| ... | ... | |
| 275 | 337 |
### delete the temporary files |
| 276 | 338 |
unlink_.gislock() |
| 277 | 339 |
system(paste("rm -frR ",tf,sep=""))
|
| 278 |
} |
|
| 279 | 340 |
|
| 280 | 341 |
|
| 281 |
########################################### |
|
| 282 |
### Define a wrapper function that will call the two functions above (gridding and QA-handling) for a single tile-date |
|
| 342 |
## Confirm that the file has the correct attributes, otherwise delete it |
|
| 343 |
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
|
| 344 |
## confirm it has all 'final variables as specified above" |
|
| 345 |
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
|
|
| 283 | 346 |
|
| 284 |
mod06<-function(date,tile){
|
|
| 285 |
print(paste("Processing date ",date," for tile",tile))
|
|
| 286 |
##################################################### |
|
| 287 |
## Run the gridding procedure |
|
| 288 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
|
| 289 |
|
|
| 290 |
##find swaths in region from sqlite database for the specified date/tile |
|
| 291 |
if(verbose) print("Accessing swath ID's from database")
|
|
| 292 |
con=dbConnect("SQLite", dbname = db)
|
|
| 293 |
fs=dbGetQuery(con,paste("SELECT * from swath_geo
|
|
| 294 |
WHERE east>=",tile_bb$lon_min," AND |
|
| 295 |
west<=",tile_bb$lon_max," AND |
|
| 296 |
north>=",tile_bb$lat_min," AND |
|
| 297 |
south<=",tile_bb$lat_max," AND |
|
| 298 |
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND |
|
| 299 |
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j")) |
|
| 300 |
)) |
|
| 301 |
con=dbDisconnect(con) |
|
| 302 |
fs$id=substr(fs$id,7,19) |
|
| 303 |
if(verbose) print(paste("###############",nrow(fs)," swath IDs recieved from database"))
|
|
| 304 |
|
|
| 305 |
## find the swaths on disk (using datadir) |
|
| 306 |
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T) |
|
| 307 |
|
|
| 308 |
## Run the gridding procedure |
|
| 309 |
lapply(swaths,swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max)) |
|
| 310 |
|
|
| 311 |
## confirm at least one file for this date is present |
|
| 312 |
outfiles=paste(tempdir(),"/",basename(swaths),sep="") |
|
| 313 |
if(!any(file.exists(outfiles))) {
|
|
| 314 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
|
| 315 |
q("no",status=0)
|
|
| 316 |
} |
|
| 317 |
|
|
| 318 |
##################################################### |
|
| 319 |
## Process the gridded files |
|
| 320 |
|
|
| 321 |
## run the mod06 processing for this date |
|
| 322 |
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
|
| 323 |
loadcloud(date,swaths=swaths,ncfile=ncfile) |
|
| 324 |
|
|
| 325 |
## Confirm that the file has the correct attributes, otherwise delete it |
|
| 326 |
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
|
| 327 |
## confirm it has all 'final variables as specified above" |
|
| 328 |
fvar=all(finalvars%in%do.call(c,strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")))
|
|
| 329 |
|
|
| 330 |
if(!ntime==1&fvar) {
|
|
| 347 |
if(ntime!=1|!fvar) {
|
|
| 331 | 348 |
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
| 332 | 349 |
file.remove(ncfile) |
| 333 | 350 |
} |
| 334 | 351 |
|
| 335 | 352 |
## print out some info |
| 336 |
print(paste(" ################################################################### Finished ",date,"
|
|
| 337 |
################################################################")) |
|
| 338 |
} |
|
| 339 |
|
|
| 340 |
## test it |
|
| 341 |
##date=notdone[1] |
|
| 342 |
mod06(date,tile) |
|
| 343 |
|
|
| 344 |
## run it for all dates - Use this if running on a workstation/server (otherwise use qsub) |
|
| 345 |
#mclapply(notdone,mod06,tile,mc.cores=ncores) # use ncores/2 because system() commands can add second process for each spawned R |
|
| 346 |
#foreach(i=notdone[1:3],.packages=(.packages())) %dopar% mod06(i,tile) |
|
| 347 |
#foreach(i=1:20) %dopar% print(i) |
|
| 353 |
print(paste(" ################################################################### Finished ",date,
|
|
| 354 |
"################################################################")) |
|
| 348 | 355 |
|
| 349 | 356 |
## delete old files |
| 350 | 357 |
system("cleartemp")
|
| 351 | 358 |
|
| 359 |
## quit |
|
| 352 | 360 |
q("no",status=0)
|
Also available in: Unified diff
restructured MOD06 processing routine to eliminate functions for more transparent (linear) processing. Also updated to run on either pleiades or litoria