Revision e5c2e69b
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD35_L2_process.r | ||
|---|---|---|
| 198 | 198 |
i=1 |
| 199 | 199 |
system(paste("gdalinfo ",outfiles[19]))
|
| 200 | 200 |
d=lapply(outfiles,function(r) raster(r)) |
| 201 |
summary(d[[5]])
|
|
| 201 |
summary(d[[6]])
|
|
| 202 | 202 |
} |
| 203 | 203 |
#system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
| 204 | 204 |
|
| ... | ... | |
| 306 | 306 |
## use r.series to find minimum |
| 307 | 307 |
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
| 308 | 308 |
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
| 309 |
|
|
| 310 | 309 |
## select cloud observation with lowest sensor zenith for day and night |
| 311 | 310 |
system( |
| 312 | 311 |
paste("r.mapcalc <<EOF
|
| ... | ... | |
| 314 | 313 |
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
| 315 | 314 |
)) |
| 316 | 315 |
|
| 316 |
execGRASS("r.null",map="CMday_daily",setnull="255") ; print("")
|
|
| 317 |
execGRASS("r.null",map="CMnight_daily",setnull="255") ; print("")
|
|
| 318 |
|
|
| 317 | 319 |
if(plot){
|
| 318 | 320 |
ps=1:nfs |
| 319 |
ps=c(12,14,17)
|
|
| 321 |
ps=c(10,11,13,14)
|
|
| 320 | 322 |
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZday_",i,sep="")))))
|
| 321 | 323 |
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMday_",i,sep="")))))
|
| 322 | 324 |
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
| ... | ... | |
| 329 | 331 |
|
| 330 | 332 |
### Write the files to a netcdf file |
| 331 | 333 |
## create image group to facilitate export as multiband netcdf |
| 332 |
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily")) ; print("")
|
|
| 334 |
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily"),flags=c("quiet")) ; print("")
|
|
| 333 | 335 |
|
| 334 | 336 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it |
| 335 |
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
|
| 337 |
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("verbose"),
|
|
| 336 | 338 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 337 | 339 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 338 | 340 |
|
| ... | ... | |
| 353 | 355 |
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 354 | 356 |
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep="")) |
| 355 | 357 |
## add other attributes |
| 358 |
## need to delete _FillValue becuase r.out.gdal incorrectly calls zero values missing if there are no other missing values in the raster. |
|
| 359 |
## so need to delete then re-add. If you just change the value, ncatted will change the values in the raster in addition to the attribute. |
|
| 356 | 360 |
system(paste(ncopath,"ncrename -v Band1,CMday -v Band2,CMnight ",ncfile,sep="")) |
| 357 | 361 |
system(paste(ncopath,"ncatted ", |
| 358 | 362 |
" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ", |
| 359 | 363 |
" -a missing_value,CMday,o,b,255 ", |
| 360 |
" -a _FillValue,CMday,o,b,255 ",
|
|
| 364 |
" -a _FillValue,CMday,d,, ",
|
|
| 361 | 365 |
" -a valid_range,CMday,o,b,\"0,3\" ", |
| 362 | 366 |
" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ", |
| 363 | 367 |
" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ", |
| 364 | 368 |
" -a missing_value,CMnight,o,b,255 ", |
| 365 |
" -a _FillValue,CMnight,o,b,255 ",
|
|
| 369 |
" -a _FillValue,CMnight,d,, ",
|
|
| 366 | 370 |
" -a valid_range,CMnight,o,b,\"0,3\" ", |
| 367 | 371 |
" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ", |
| 368 | 372 |
ncfile,sep="")) |
| 369 |
#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep="")) |
|
| 373 |
## add the fillvalue attribute back (without changing the actual values) |
|
| 374 |
system(paste(ncopath,"ncatted -a _FillValue,CMday,o,b,255 ",ncfile,sep="")) |
|
| 375 |
system(paste(ncopath,"ncatted -a _FillValue,CMnight,o,b,255 ",ncfile,sep="")) |
|
| 370 | 376 |
|
| 371 | 377 |
|
| 372 | 378 |
## Confirm that the file has the correct attributes, otherwise delete it |
| climate/procedures/Pleiades_MOD35.R | ||
|---|---|---|
| 103 | 103 |
### report on what has already been processed |
| 104 | 104 |
print(paste(sum(!proclist$done)," out of ",nrow(proclist)," (",round(100*sum(!proclist$done)/nrow(proclist),2),"%) remain"))
|
| 105 | 105 |
stem(table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done])) |
| 106 |
table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done]) |
|
| 106 |
#table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done])
|
|
| 107 | 107 |
table(table(tile=proclist$tile[!proclist$done],year=proclist$year[!proclist$done])) |
| 108 | 108 |
|
| 109 | 109 |
### explore tile counts |
| ... | ... | |
| 127 | 127 |
i=2 |
| 128 | 128 |
time1=system.time(system(paste("Rscript --verbose ",script," --date ",proclist$date[i]," --verbose T --tile ",proclist$tile[i],sep="")))
|
| 129 | 129 |
hours=round(length(proclist$date[tp])*142/60/60) |
| 130 |
hours=round(length(proclist$date[tp])*time1[3]/60/60,1) |
|
| 131 |
hours/240
|
|
| 130 |
hours=round(length(proclist$date[tp])*time1[3]/60/60,1); hours
|
|
| 131 |
hours/400
|
|
| 132 | 132 |
print(paste("Based on runtime of previous command, it will take",hours," hours to process the full set"))
|
| 133 | 133 |
} |
| 134 | 134 |
|
| ... | ... | |
| 136 | 136 |
### qsub script |
| 137 | 137 |
cat(paste("
|
| 138 | 138 |
#PBS -S /bin/bash |
| 139 |
#PBS -l select=28:ncpus=8:mpiprocs=8
|
|
| 139 |
#PBS -l select=50:ncpus=8:mpiprocs=8
|
|
| 140 | 140 |
##PBS -l select=100:ncpus=8:mpiprocs=8 |
| 141 | 141 |
##PBS -l walltime=8:00:00 |
| 142 | 142 |
#PBS -l walltime=2:00:00 |
| ... | ... | |
| 148 | 148 |
#PBS -V |
| 149 | 149 |
|
| 150 | 150 |
#CORES=800 |
| 151 |
CORES=224
|
|
| 151 |
CORES=400
|
|
| 152 | 152 |
|
| 153 | 153 |
HDIR=/u/armichae/pr/ |
| 154 | 154 |
source $HDIR/etc/environ.sh |
| ... | ... | |
| 179 | 179 |
### Now submit the script to generate the climatologies |
| 180 | 180 |
|
| 181 | 181 |
## report 'mostly' finished tiles |
| 182 |
## this relyies on proclist above so be sure to update above before running
|
|
| 182 |
## this relies on proclist above so be sure to update above before running |
|
| 183 | 183 |
md=table(tile=proclist$tile[!proclist$done],year=proclist$year[!proclist$done]) |
| 184 | 184 |
mdt=names(md[md<10,]) |
| 185 | 185 |
tiles=mdt |
| ... | ... | |
| 204 | 204 |
file=paste("notdone_climate.txt",sep=""),row.names=F,col.names=F,quote=F)
|
| 205 | 205 |
|
| 206 | 206 |
## delay start until previous jobs have finished? |
| 207 |
delay=F
|
|
| 207 |
delay=T
|
|
| 208 | 208 |
## check running jobs to get JobID of job you want to wait for |
| 209 |
system("qstat -u awilso10")
|
|
| 209 |
system("qstat -u awilso10",intern=T)
|
|
| 210 | 210 |
## enter JobID here: |
| 211 |
job="881394.pbspl1.nas.nasa.gov"
|
|
| 211 |
job="2031668.pbspl1.nas.nasa.gov"
|
|
| 212 | 212 |
|
| 213 | 213 |
### qsub script |
| 214 | 214 |
cat(paste("
|
| 215 | 215 |
#PBS -S /bin/bash |
| 216 |
#PBS -l select=10:ncpus=8:mem=94
|
|
| 216 |
#PBS -l select=4:ncpus=8:mem=94
|
|
| 217 | 217 |
#PBS -l walltime=2:00:00 |
| 218 | 218 |
#PBS -j n |
| 219 | 219 |
#PBS -m be |
| ... | ... | |
| 224 | 224 |
#PBS -V |
| 225 | 225 |
",if(delay) paste("#PBS -W depend=afterany:",job,sep=""),"
|
| 226 | 226 |
|
| 227 |
CORES=80
|
|
| 227 |
CORES=32
|
|
| 228 | 228 |
HDIR=/u/armichae/pr/ |
| 229 | 229 |
source $HDIR/etc/environ.sh |
| 230 | 230 |
source /pleiades/u/awilso10/environ.sh |
Also available in: Unified diff
found and fixed bug in r.out.gdal export of netcdf files with no missing data. The missing data value was set to 0 instead of the 255 that was specified. This resulted in all 0s to be removed from these day-tiles. Fixed by deleting and resetting the attribute _FillValue using ncatted.