Revision f27d277d
Added by Adam Wilson about 11 years ago
| 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="20010410"
|
|
| 38 |
date="20000410"
|
|
| 39 | 39 |
tile="h11v08" |
| 40 | 40 |
platform="pleiades" |
| 41 | 41 |
verbose=T |
| ... | ... | |
| 68 | 68 |
|
| 69 | 69 |
if(platform=="litoria"){ #if running on local server, use different paths
|
| 70 | 70 |
## specify working directory |
| 71 |
setwd("~/acrobates/projects/interp")
|
|
| 71 |
setwd("~/acrobates/adamw/projects/interp")
|
|
| 72 | 72 |
gisBase="/usr/lib/grass64" |
| 73 | 73 |
## location of MOD06 files |
| 74 |
datadir="~/acrobates/projects/interp/data/modis/mod35" |
|
| 74 |
datadir="~/acrobates/adamw/projects/interp/data/modis/mod35"
|
|
| 75 | 75 |
## path to some executables |
| 76 | 76 |
ncopath="" |
| 77 | 77 |
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif" |
| 78 | 78 |
## path to swath database |
| 79 |
db="/home/adamw/acrobates/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
|
|
| 79 |
db="~/acrobates/adamw/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
|
|
| 80 | 80 |
## get grid file |
| 81 |
td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
|
| 81 |
td=raster(paste("~/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
|
| 82 | 82 |
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 " |
| 83 | 83 |
} |
| 84 | 84 |
|
| ... | ... | |
| 100 | 100 |
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
| 101 | 101 |
|
| 102 | 102 |
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end. |
| 103 |
finalvars=c("CLD","CLD2")
|
|
| 103 |
finalvars=c("PClear")
|
|
| 104 | 104 |
|
| 105 | 105 |
|
| 106 | 106 |
##################################################### |
| 107 | 107 |
##find swaths in region from sqlite database for the specified date/tile |
| 108 | 108 |
if(verbose) print("Accessing swath ID's from database")
|
| 109 | 109 |
con=dbConnect("SQLite", dbname = db)
|
| 110 |
fs=dbGetQuery(con,paste("SELECT * from swath_geo
|
|
| 110 |
fs=dbGetQuery(con,paste("SELECT * from swath_geo6
|
|
| 111 | 111 |
WHERE east>=",tile_bb$lon_min," AND |
| 112 | 112 |
west<=",tile_bb$lon_max," AND |
| 113 | 113 |
north>=",tile_bb$lat_min," AND |
| ... | ... | |
| 209 | 209 |
|
| 210 | 210 |
## Define region by importing one MOD11A1 raster. |
| 211 | 211 |
print("Import one MOD11A1 raster to define grid")
|
| 212 |
if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 213 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 214 |
if(platform=="litoria") |
|
| 215 |
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("overwrite","o"))
|
|
| 216 |
|
|
| 217 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 218 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 219 |
|
|
| 212 |
if(platform=="pleiades") {
|
|
| 213 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 214 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 215 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 216 |
} |
|
| 217 |
|
|
| 218 |
if(platform=="litoria"){
|
|
| 219 |
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
|
| 220 |
output="modisgrid",flags=c("overwrite","o"))
|
|
| 221 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 222 |
} |
|
| 220 | 223 |
## Identify which files to process |
| 221 | 224 |
tfs=basename(swaths) |
| 222 | 225 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
| ... | ... | |
| 231 | 234 |
## Cloud Mask |
| 232 | 235 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
|
| 233 | 236 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 234 |
## extract cloudy and 'probably/confidently clear' pixels |
|
| 237 |
## QA ## extract first bit to keep only "useful" values of cloud mask |
|
| 238 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
|
|
| 239 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 240 |
## extract first two bits of cloud mask |
|
| 235 | 241 |
system(paste("r.mapcalc <<EOF
|
| 236 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
|
| 237 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) >= 2 |
|
| 242 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0) |
|
| 243 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999) |
|
| 244 |
Pclear1_",i," = if(CM_cloud_",i,"==0,0,if(CM_cloud_",i,"==1,66,if(CM_cloud_",i,"==2,95,if(CM_cloud_",i,"==3,99,-9999)))) |
|
| 245 |
Pclear_",i," = if(QA_useful_",i,"==1,Pclear1_",i,",-9999) |
|
| 238 | 246 |
CM_path_",i," = ((CM1_",i," / 2^6) % 2^2) |
| 239 |
CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2) |
|
| 240 | 247 |
EOF",sep="")) |
| 241 | 248 |
|
| 242 |
## QA |
|
| 243 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
|
|
| 244 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 245 |
## QA_CER |
|
| 246 |
# system(paste("r.mapcalc <<EOF
|
|
| 247 |
# QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1 |
|
| 248 |
# QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )>=2 |
|
| 249 |
# QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0 |
|
| 250 |
# QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1 |
|
| 251 |
# QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )>=2 |
|
| 252 |
#EOF",sep="")) |
|
| 253 |
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1 |
|
| 254 |
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3 |
|
| 255 |
|
|
| 256 |
|
|
| 249 |
execGRASS("r.null",map=paste("Pclear_",i,sep=""),setnull="-9999")
|
|
| 250 |
execGRASS("r.null",map=paste("CM_cloud_",i,sep=""),setnull="-9999")
|
|
| 251 |
|
|
| 252 |
|
|
| 257 | 253 |
} #end loop through sub daily files |
| 258 | 254 |
|
| 259 | 255 |
#### Now generate daily averages (or maximum in case of cloud flag) |
| 260 | 256 |
|
| 261 | 257 |
system(paste("r.mapcalc <<EOF
|
| 262 |
CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),-9999,CM_cloud_",1:nfs,")",sep="",collapse=","),")))
|
|
| 263 |
CLD2_daily=int((min(",paste("if(isnull(CM_cloud2_",1:nfs,"),-9999,CM_cloud2_",1:nfs,")",sep="",collapse=","),")))
|
|
| 258 |
Pclear_daily=int((min(",paste("if(isnull(Pclear_",1:nfs,"),9999,Pclear_",1:nfs,")",sep="",collapse=","),")))
|
|
| 264 | 259 |
EOF",sep="")) |
| 265 | 260 |
|
| 266 |
execGRASS("r.null",map="CLD_daily",setnull="-9999")
|
|
| 267 |
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
|
|
| 261 |
# CLD_daily=int((min(",paste("if(isnull(CM_cloud_",1:nfs,"),9999,CM_cloud_",1:nfs,")",sep="",collapse=","),")))
|
|
| 262 |
|
|
| 263 |
## reset null values |
|
| 264 |
#execGRASS("r.null",map="CLD_daily",setnull="9999")
|
|
| 265 |
execGRASS("r.null",map="Pclear_daily",setnull="9999")
|
|
| 268 | 266 |
|
| 269 | 267 |
|
| 270 | 268 |
### Write the files to a netcdf file |
| 271 | 269 |
## create image group to facilitate export as multiband netcdf |
| 272 |
execGRASS("i.group",group="mod06",input=c("CLD_daily","CLD2_daily")) ; print("")
|
|
| 270 |
execGRASS("i.group",group="mod35",input=c("Pclear_daily")) ; print("")
|
|
| 273 | 271 |
|
| 274 | 272 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it |
| 275 |
execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),
|
|
| 273 |
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
|
| 276 | 274 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 277 | 275 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 278 | 276 |
|
| ... | ... | |
| 293 | 291 |
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 294 | 292 |
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep="")) |
| 295 | 293 |
## add other attributes |
| 296 |
system(paste(ncopath,"ncrename -v Band1,CLD -v Band2,CLD2 ",ncfile,sep="")) |
|
| 297 |
system(paste(ncopath,"ncatted -a scale_factor,CLD,o,d,1 -a units,CLD,o,c,\"none\" -a missing_value,CLD,o,d,-32768 -a long_name,CLD,o,c,\"Cloud Mask\" ",ncfile,sep="")) |
|
| 298 |
system(paste(ncopath,"ncatted -a scale_factor,CLD2,o,d,1 -a units,CLD2,o,c,\"none\" -a missing_value,CLD2,o,d,-32768 -a long_name,CLD2,o,c,\"Cloud Mask Flag\" ",ncfile,sep="")) |
|
| 294 |
system(paste(ncopath,"ncrename -v Band1,PClear ",ncfile,sep="")) |
|
| 295 |
system(paste(ncopath,"ncatted -a scale_factor,PClear,o,d,1 -a units,PClear,o,c,\"Probability (%)\" -a missing_value,PClear,o,d,255 -a long_name,PClear,o,c,\"Probability of Clear Sky\" ",ncfile,sep="")) |
|
| 299 | 296 |
|
| 300 | 297 |
# system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep="")) |
| 301 | 298 |
|
| ... | ... | |
| 320 | 317 |
"################################################################")) |
| 321 | 318 |
|
| 322 | 319 |
## delete old files |
| 323 |
system("cleartemp")
|
|
| 320 |
#system("cleartemp")
|
|
| 324 | 321 |
|
| 325 | 322 |
## quit |
| 326 | 323 |
q("no",status=0)
|
Also available in: Unified diff
Simplified MOD35 script and converted output to single band p(Clear) rather than multiple bands