Revision 3bb88310
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD09_cloud.r | ||
|---|---|---|
| 1 | 1 |
################################################################################### |
| 2 |
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
|
| 2 |
### R code to aquire and process MOD35_L2 cloud data from the MODIS platform
|
|
| 3 | 3 |
|
| 4 | 4 |
|
| 5 | 5 |
# Redirect all warnings to stderr() |
| ... | ... | |
| 17 | 17 |
require(raster) |
| 18 | 18 |
require(rgdal) |
| 19 | 19 |
require(spgrass6) |
| 20 |
require(RSQLite) |
|
| 21 | 20 |
|
| 22 | 21 |
## get options |
| 23 | 22 |
opta <- getopt(matrix(c( |
| ... | ... | |
| 33 | 32 |
q(status=1); |
| 34 | 33 |
} |
| 35 | 34 |
|
| 35 |
testing=T |
|
| 36 |
platform="pleiades" |
|
| 36 | 37 |
|
| 37 | 38 |
## default date and tile to play with (will be overwritten below when running in batch) |
| 38 |
date="20030102" |
|
| 39 |
tile="h11v08" |
|
| 40 |
platform="pleiades" |
|
| 41 |
verbose=T |
|
| 39 |
if(testing){
|
|
| 40 |
date="20090129" |
|
| 41 |
tile="h11v08" |
|
| 42 |
tile="h17v00" |
|
| 43 |
verbose=T |
|
| 44 |
} |
|
| 42 | 45 |
|
| 43 | 46 |
## now update using options if given |
| 44 |
date=opta$date |
|
| 45 |
tile=opta$tile |
|
| 46 |
verbose=opta$verbose #print out extensive information for debugging? |
|
| 47 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
|
| 48 |
|
|
| 49 |
|
|
| 47 |
if(!testing){
|
|
| 48 |
date=opta$date |
|
| 49 |
tile=opta$tile |
|
| 50 |
verbose=opta$verbose #print out extensive information for debugging? |
|
| 51 |
} |
|
| 50 | 52 |
## get year and doy from date |
| 51 | 53 |
year=format(as.Date(date,"%Y%m%d"),"%Y") |
| 52 | 54 |
doy=format(as.Date(date,"%Y%m%d"),"%j") |
| 53 | 55 |
|
| 54 | 56 |
if(platform=="pleiades"){
|
| 55 |
## location of MOD06 files
|
|
| 56 |
datadir=paste("/nobackupp4/datapool/modis/MOD06_L2.005/",year,"/",doy,"/",sep="")
|
|
| 57 |
## location of MOD35 files
|
|
| 58 |
datadir=paste("/nobackupp4/datapool/modis/MOD09GA.005/",year,"/",doy,"/",sep="")
|
|
| 57 | 59 |
## path to some executables |
| 58 | 60 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
| 59 |
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif" |
|
| 60 |
## path to swath database |
|
| 61 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 62 | 61 |
## specify working directory |
| 63 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 62 |
outdir=paste("/nobackupp1/awilso10/mod09/daily/",tile,"/",sep="") #directory for separate daily files
|
|
| 63 |
basedir="/nobackupp1/awilso10/mod09/" #directory to hold files temporarily before transferring to lou |
|
| 64 |
setwd(tempdir()) |
|
| 65 |
## grass database |
|
| 64 | 66 |
gisBase="/u/armichae/pr/grass-6.4.2/" |
| 65 | 67 |
## path to MOD11A1 file for this tile to align grid/extent |
| 66 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
|
| 67 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
|
| 68 |
gridfile=list.files("/nobackupp1/awilso10/mod35/MODTILES/",pattern=tile,full=T)[1]
|
|
| 69 |
td=raster(gridfile)
|
|
| 68 | 70 |
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 " |
| 69 | 71 |
} |
| 70 | 72 |
|
| 71 |
if(platform=="litoria"){ #if running on local server, use different paths
|
|
| 72 |
## specify working directory |
|
| 73 |
setwd("~/acrobates/projects/interp")
|
|
| 74 |
gisBase="/usr/lib/grass64" |
|
| 75 |
## location of MOD06 files |
|
| 76 |
datadir="~/acrobates/projects/interp/data/modis/Venezuela/MOD06" |
|
| 77 |
## path to some executables |
|
| 78 |
ncopath="" |
|
| 79 |
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif" |
|
| 80 |
## path to swath database |
|
| 81 |
db="/home/adamw/acrobates/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db" |
|
| 82 |
## get grid file |
|
| 83 |
td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
|
| 84 |
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 " |
|
| 85 |
} |
|
| 86 |
|
|
| 87 |
|
|
| 88 | 73 |
### print some status messages |
| 89 | 74 |
if(verbose) print(paste("Processing tile",tile," for date",date))
|
| 90 | 75 |
|
| 91 | 76 |
## load tile information and get bounding box |
| 92 |
load(file="modlandTiles.Rdata") |
|
| 77 |
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata")
|
|
| 93 | 78 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
| 94 |
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner |
|
| 95 |
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner |
|
| 96 |
|
|
| 97 |
|
|
| 98 |
## vars to process |
|
| 99 |
vars=as.data.frame(matrix(c( |
|
| 100 |
"Cloud_Effective_Radius", "CER", |
|
| 101 |
"Cloud_Effective_Radius_Uncertainty", "CERU", |
|
| 102 |
"Cloud_Optical_Thickness", "COT", |
|
| 103 |
"Cloud_Optical_Thickness_Uncertainty", "COTU", |
|
| 104 |
# "Cloud_Water_Path", "CWP", |
|
| 105 |
# "Cloud_Water_Path_Uncertainty", "CWPU", |
|
| 106 |
# "Cloud_Phase_Optical_Properties", "CPOP", |
|
| 107 |
# "Cloud_Multi_Layer_Flag", "CMLF", |
|
| 108 |
"Cloud_Mask_1km", "CM1", |
|
| 109 |
"Quality_Assurance_1km", "QA"), |
|
| 110 |
byrow=T,ncol=2,dimnames=list(1:6,c("variable","varid"))),stringsAsFactors=F)
|
|
| 111 |
|
|
| 112 |
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end. |
|
| 113 |
finalvars=c("CER","COT","CLD")
|
|
| 114 |
|
|
| 115 | 79 |
|
| 116 | 80 |
##################################################### |
| 117 |
##find swaths in region from sqlite database for the specified date/tile |
|
| 118 |
if(verbose) print("Accessing swath ID's from database")
|
|
| 119 |
con=dbConnect("SQLite", dbname = db)
|
|
| 120 |
fs=dbGetQuery(con,paste("SELECT * from swath_geo
|
|
| 121 |
WHERE east>=",tile_bb$lon_min," AND |
|
| 122 |
west<=",tile_bb$lon_max," AND |
|
| 123 |
north>=",tile_bb$lat_min," AND |
|
| 124 |
south<=",tile_bb$lat_max," AND |
|
| 125 |
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND |
|
| 126 |
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j")) |
|
| 127 |
)) |
|
| 128 |
con=dbDisconnect(con) |
|
| 129 |
fs$id=substr(fs$id,7,19) |
|
| 130 |
## find the swaths on disk (using datadir) |
|
| 131 |
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T) |
|
| 132 |
|
|
| 133 |
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk")) |
|
| 134 |
|
|
| 135 |
|
|
| 136 |
############################################################################ |
|
| 137 |
############################################################################ |
|
| 138 |
### Use the HEG tool to grid all available swath data for this date-tile |
|
| 139 |
for(file in swaths){
|
|
| 140 |
## Function to generate hegtool parameter file for multi-band HDF-EOS file |
|
| 141 |
print(paste("Starting file",basename(file)))
|
|
| 142 |
outfile=paste(tempdir(),"/",basename(file),sep="") |
|
| 143 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 144 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 145 |
grp=paste("
|
|
| 146 |
BEGIN |
|
| 147 |
INPUT_FILENAME=",file," |
|
| 148 |
OBJECT_NAME=mod06 |
|
| 149 |
FIELD_NAME=",vars$variable,"| |
|
| 150 |
BAND_NUMBER = 1 |
|
| 151 |
OUTPUT_PIXEL_SIZE_X=1000 |
|
| 152 |
OUTPUT_PIXEL_SIZE_Y=1000 |
|
| 153 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," ) |
|
| 154 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," ) |
|
| 155 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
|
| 156 |
RESAMPLING_TYPE =NN |
|
| 157 |
OUTPUT_PROJECTION_TYPE = SIN |
|
| 158 |
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ) |
|
| 159 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
|
| 160 |
ELLIPSOID_CODE = WGS84 |
|
| 161 |
OUTPUT_TYPE = HDFEOS |
|
| 162 |
OUTPUT_FILENAME = ",outfile," |
|
| 163 |
END |
|
| 164 |
|
|
| 165 |
",sep="") |
|
| 166 |
|
|
| 167 |
## if any remnants from previous runs remain, delete them |
|
| 168 |
if(length(list.files(tempdir(),pattern=basename(file)))>0) |
|
| 169 |
file.remove(list.files(tempdir(),pattern=basename(file),full=T)) |
|
| 170 |
## write it to a file |
|
| 171 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
|
| 172 |
## now run the swath2grid tool |
|
| 173 |
## write the gridded file |
|
| 174 |
log=system(paste("(",swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
|
| 175 |
## clean up temporary files in working directory |
|
| 176 |
# file.remove(list.files(pattern= |
|
| 177 |
# paste("filetable.temp_",
|
|
| 178 |
# as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process |
|
| 179 |
if(verbose) print(log) |
|
| 180 |
print(paste("Finished gridding ", file))
|
|
| 181 |
} #end looping over swaths |
|
| 182 |
|
|
| 183 |
######################## |
|
| 184 |
## confirm at least one file for this date is present. If not, quit. |
|
| 185 |
outfiles=paste(tempdir(),"/",basename(swaths),sep="") |
|
| 186 |
if(!any(file.exists(outfiles))) {
|
|
| 187 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
|
| 188 |
q("no",status=0)
|
|
| 189 |
} |
|
| 190 |
|
|
| 191 |
##################################################### |
|
| 192 |
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths |
|
| 193 | 81 |
|
| 194 |
## Identify output file |
|
| 195 |
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
|
| 196 |
|
|
| 197 | 82 |
## function to convert binary to decimal to assist in identifying correct values |
| 198 | 83 |
## this is helpful when defining QA handling below, but isn't used in processing |
| 199 | 84 |
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html |
| ... | ... | |
| 210 | 95 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar |
| 211 | 96 |
if(!file.exists(tf)) dir.create(tf) |
| 212 | 97 |
## create output directory if needed |
| 98 |
## Identify output file |
|
| 99 |
ncfile=paste(outdir,"MOD09cloud_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file |
|
| 213 | 100 |
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T) |
| 214 |
|
|
| 101 |
|
|
| 215 | 102 |
## set up temporary grass instance for this PID |
| 216 | 103 |
if(verbose) print(paste("Set up temporary grass session in",tf))
|
| 217 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 104 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod09",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 218 | 105 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 219 | 106 |
|
| 220 | 107 |
## Define region by importing one MOD11A1 raster. |
| 221 | 108 |
print("Import one MOD11A1 raster to define grid")
|
| 222 |
if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 223 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 224 |
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("overwrite","o"))
|
|
| 225 |
|
|
| 226 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 227 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 109 |
if(platform=="pleiades") {
|
|
| 110 |
execGRASS("r.in.gdal",input=td@file@name,output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 111 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 112 |
} |
|
| 113 |
|
|
| 114 |
if(platform=="litoria"){
|
|
| 115 |
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
|
| 116 |
output="modisgrid",flags=c("overwrite","o"))
|
|
| 117 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
| 118 |
} |
|
| 228 | 119 |
|
| 229 | 120 |
## Identify which files to process |
| 230 |
tfs=basename(swaths) |
|
| 231 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
|
| 232 |
tfs=tfs[tfs%in%list.files(tempdir())] |
|
| 121 |
tfs=unique(sub("CM_|QA_|SenZen_|SolZen_","",basename(outfiles)))
|
|
| 122 |
#tfs=list.files(tempdir(),pattern="temp.*hdf") |
|
| 233 | 123 |
nfs=length(tfs) |
| 234 | 124 |
if(verbose) print(paste(nfs,"swaths available for processing")) |
| 235 | 125 |
|
| 236 | 126 |
## loop through scenes and process QA flags |
| 237 | 127 |
for(i in 1:nfs){
|
| 238 |
file=paste(tempdir(),"/",tfs[i],sep="") |
|
| 128 |
bfile=tfs[i] |
|
| 129 |
## Read in the data from the HDFs |
|
| 239 | 130 |
## Cloud Mask |
| 240 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
| 131 |
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
|
| 241 | 132 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 242 |
## extract cloudy and 'probably/confidently clear' pixels |
|
| 243 |
system(paste("r.mapcalc <<EOF
|
|
| 244 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
|
| 245 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) > 2 |
|
| 246 |
CM_path_",i," = ((CM1_",i," / 2^6) % 2^2) |
|
| 247 |
CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2) |
|
| 248 |
EOF",sep="")) |
|
| 249 |
## Set CM_cloud2 to null if it is "01" (uncertain) |
|
| 250 |
# execGRASS("r.null",map=paste("CM_cloud2_",i,sep=""),setnull="-9999,1")
|
|
| 251 |
|
|
| 252 |
## QA |
|
| 253 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
|
| 133 |
## QA ## extract first bit to keep only "useful" values of cloud mask |
|
| 134 |
execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),
|
|
| 254 | 135 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 255 |
## QA_CER |
|
| 256 |
system(paste("r.mapcalc <<EOF
|
|
| 257 |
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1 |
|
| 258 |
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )>=2 |
|
| 259 |
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0 |
|
| 260 |
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1 |
|
| 261 |
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )>=2 |
|
| 262 |
EOF",sep="")) |
|
| 263 |
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1 |
|
| 264 |
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3 |
|
| 265 |
|
|
| 266 |
## Optical Thickness |
|
| 267 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
|
| 268 |
output=paste("COT_",i,sep=""),
|
|
| 269 |
title="cloud_effective_radius", |
|
| 270 |
flags=c("overwrite","o")) ; print("")
|
|
| 271 |
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
|
| 272 |
## keep only positive COT values where quality is 'useful' and '>= good' & scale to real units |
|
| 273 |
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,",null())\"",sep=""))
|
|
| 274 |
## set null COT to 0 in clear-sky pixels |
|
| 275 |
system(paste("r.mapcalc \"COT2_",i,"=if(isnull(COT_",i,")&&CM_clear_",i,"==1,0,COT_",i,")\"",sep=""))
|
|
| 136 |
## Sensor Zenith ## extract first bit to keep only "low angle" observations |
|
| 137 |
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
|
| 138 |
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 276 | 139 |
|
| 277 |
## Effective radius ## |
|
| 278 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
|
| 279 |
output=paste("CER_",i,sep=""),
|
|
| 280 |
title="cloud_effective_radius", |
|
| 281 |
flags=c("overwrite","o")) ; print("")
|
|
| 282 |
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
|
| 283 |
## keep only positive CER values where quality is 'useful' and '>= good' & scale to real units |
|
| 284 |
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,",null())\"",sep=""))
|
|
| 285 |
## set null CER to 0 in clear-sky pixels |
|
| 286 |
system(paste("r.mapcalc \"CER2_",i,"=if(isnull(CER_",i,")&&CM_clear_",i,"==1,0,CER_",i,")\"",sep=""))
|
|
| 287 |
|
|
| 288 |
## Cloud Water Path |
|
| 289 |
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
|
| 290 |
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
|
| 291 |
# flags=c("overwrite","o")) ; print("")
|
|
| 292 |
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
|
| 293 |
# system(paste("r.mapcalc \"CWP_",i,"=if(CWP_",i,"<0,null(),CWP_",i,")\"",sep=""))
|
|
| 294 |
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units |
|
| 295 |
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,",null())\"",sep=""))
|
|
| 296 |
## set CER to 0 in clear-sky pixels |
|
| 297 |
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
|
| 140 |
## check for interpolation artefacts |
|
| 141 |
# execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))
|
|
| 142 |
# execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))
|
|
| 143 |
# execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))
|
|
| 144 |
|
|
| 145 |
## write out the table of weights for use in the neighborhood analysis to identify bad pixels from swtif |
|
| 146 |
p=75 #must be odd |
|
| 147 |
mat=matrix(rep(0,p*p),nrow=p) |
|
| 148 |
mat[0.5+p/2,]=1 |
|
| 149 |
cat(mat,file="weights.txt") |
|
| 150 |
execGRASS("r.neighbors",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"clump",sep=""),method="range",size=p,weight="weights.txt") # too slow!
|
|
| 151 |
system(paste("r.mapcalc \"SZ_",i,"_clump2=SZ_",i,"clump==0\"",sep=""))
|
|
| 152 |
|
|
| 153 |
# p=-50:50 |
|
| 154 |
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
| 155 |
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",min(p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
| 156 |
# vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))
|
|
| 157 |
# colnames(vals)=c("value","count")
|
|
| 158 |
# vals$count=as.numeric(as.character(vals$count)) |
|
| 159 |
# vals$value=as.numeric(as.character(vals$value)) |
|
| 160 |
# vals=na.omit(vals) |
|
| 161 |
# vals$count[vals$value==1&vals$count>10] |
|
| 162 |
# |
|
| 163 |
#plot(p~value,data=vals) |
|
| 164 |
# print(sum(vals$p[vals$p>.1])) |
|
| 298 | 165 |
|
| 299 |
} #end loop through sub daily files |
|
| 300 |
|
|
| 301 |
#### Now generate daily averages (or maximum in case of cloud flag) |
|
| 302 |
|
|
| 303 |
system(paste("r.mapcalc <<EOF
|
|
| 304 |
COT_numer=",paste("if(isnull(COT_",1:nfs,"),0,COT_",1:nfs,")",sep="",collapse="+"),"
|
|
| 305 |
COT_denom=",paste("!isnull(COT_",1:nfs,")",sep="",collapse="+"),"
|
|
| 306 |
COT_daily=int(COT_numer/COT_denom) |
|
| 307 |
COT2_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 308 |
COT2_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 309 |
COT2_daily=int(COT2_numer/COT2_denom) |
|
| 310 |
CER_numer=",paste("if(isnull(CER_",1:nfs,"),0,CER_",1:nfs,")",sep="",collapse="+"),"
|
|
| 311 |
CER_denom=",paste("!isnull(CER_",1:nfs,")",sep="",collapse="+"),"
|
|
| 312 |
CER_daily=int(CER_numer/CER_denom) |
|
| 313 |
CER2_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 314 |
CER2_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 315 |
CER2_daily=int(CER2_numer/CER2_denom) |
|
| 316 |
CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),-9999,CM_cloud_",1:nfs,")",sep="",collapse=","),")))
|
|
| 317 |
CLD2_daily=int((min(",paste("if(isnull(CM_cloud2_",1:nfs,"),-9999,CM_cloud2_",1:nfs,")",sep="",collapse=","),")))
|
|
| 166 |
## Solar Zenith ## extract first bit to keep only "low angle" observations |
|
| 167 |
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
|
| 168 |
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 169 |
## produce the summaries |
|
| 170 |
system(paste("r.mapcalc <<EOF
|
|
| 171 |
CM_fill_",i," = if(isnull(CM1_",i,"),1,0) |
|
| 172 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0) |
|
| 173 |
SZ_low_",i," = if(SZ_",i,"_clump2==0&SZ_",i,"<6000,1,0) |
|
| 174 |
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0) |
|
| 175 |
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0) |
|
| 176 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null()) |
|
| 177 |
SZday_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==1,SZ_",i,",null()) |
|
| 178 |
SZnight_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==0,SZ_",i,",null()) |
|
| 179 |
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null()) |
|
| 180 |
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null()) |
|
| 318 | 181 |
EOF",sep="")) |
| 319 | 182 |
|
| 320 |
execGRASS("r.null",map="CLD_daily",setnull="-9999")
|
|
| 321 |
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
|
|
| 183 |
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0) |
|
| 184 |
# CM_dscore_",i," = if((CM_dayflag_",i,"==0|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,if(SoZ_",i,">=8500,3,4)))) |
|
| 185 |
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4))) |
|
| 186 |
|
|
| 187 |
drawplot=F |
|
| 188 |
if(drawplot){
|
|
| 189 |
d2=stack( |
|
| 190 |
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
|
| 191 |
raster(readRAST6(paste("CM1_",i,sep=""))),
|
|
| 192 |
# raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
| 193 |
# raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
| 194 |
# raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
| 195 |
# raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
| 196 |
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
|
| 197 |
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
|
| 198 |
raster(readRAST6(paste("SZ_",i,sep=""))),
|
|
| 199 |
raster(readRAST6(paste("SZ_",i,"_clump",sep=""))),
|
|
| 200 |
raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))
|
|
| 201 |
) |
|
| 202 |
plot(d2,add=F) |
|
| 203 |
} |
|
| 204 |
|
|
| 205 |
|
|
| 206 |
} #end loop through sub daily files |
|
| 207 |
|
|
| 208 |
## select lowest view angle |
|
| 209 |
## use r.series to find minimum |
|
| 210 |
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
|
| 211 |
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
|
| 212 |
|
|
| 213 |
## select cloud observation with lowest sensor zenith for day and night |
|
| 214 |
system( |
|
| 215 |
paste("r.mapcalc <<EOF
|
|
| 216 |
CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
|
| 217 |
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
|
| 218 |
)) |
|
| 219 |
|
|
| 220 |
if(plot){
|
|
| 221 |
ps=1:nfs |
|
| 222 |
ps=c(12,14,17) |
|
| 223 |
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))
|
|
| 224 |
sz_clump=brick(lapply(ps,function(i) raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))))
|
|
| 225 |
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
| 226 |
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
|
| 227 |
library(rasterVis) |
|
| 228 |
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100)) |
|
| 229 |
levelplot(sz_clump) |
|
| 230 |
levelplot(d) |
|
| 231 |
levelplot(d2) |
|
| 232 |
} |
|
| 322 | 233 |
|
| 323 | 234 |
|
| 324 | 235 |
### Write the files to a netcdf file |
| 325 | 236 |
## create image group to facilitate export as multiband netcdf |
| 326 |
execGRASS("i.group",group="mod06",input=c("CER_daily","CER2_daily","COT_daily","COT2_daily","CLD_daily","CLD2_daily")) ; print("")
|
|
| 327 |
|
|
| 328 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
|
| 329 |
execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),
|
|
| 237 |
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily")) ; print("")
|
|
| 238 |
|
|
| 239 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it |
|
| 240 |
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
|
| 330 | 241 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 331 | 242 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 332 | 243 |
|
| ... | ... | |
| 347 | 258 |
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 348 | 259 |
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep="")) |
| 349 | 260 |
## add other attributes |
| 350 |
system(paste(ncopath,"ncrename -v Band1,CER -v Band2,CER2 -v Band3,COT -v Band4,COT2 -v Band5,CLD -v Band6,CLD2 ",ncfile,sep="")) |
|
| 351 |
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="")) |
|
| 352 |
system(paste(ncopath,"ncatted -a scale_factor,CER2,o,d,0.01 -a units,CER2,o,c,\"micron\" -a missing_value,CER2,o,d,-32768 -a long_name,CER2,o,c,\"Cloud Particle Effective Radius with clear sky set to zero\" ",ncfile,sep="")) |
|
| 353 |
|
|
| 354 |
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="")) |
|
| 355 |
system(paste(ncopath,"ncatted -a scale_factor,COT2,o,d,0.01 -a units,COT2,o,c,\"none\" -a missing_value,COT2,o,d,-32768 -a long_name,COT2,o,c,\"Cloud Optical Thickness with clear sky set to zero\" ",ncfile,sep="")) |
|
| 356 |
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="")) |
|
| 357 |
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="")) |
|
| 358 |
|
|
| 359 |
# system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep="")) |
|
| 261 |
system(paste(ncopath,"ncrename -v Band1,CMday -v Band2,CMnight ",ncfile,sep="")) |
|
| 262 |
system(paste(ncopath,"ncatted ", |
|
| 263 |
" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ", |
|
| 264 |
" -a missing_value,CMday,o,b,255 ", |
|
| 265 |
" -a _FillValue,CMday,o,b,255 ", |
|
| 266 |
" -a valid_range,CMday,o,b,\"0,3\" ", |
|
| 267 |
" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ", |
|
| 268 |
" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ", |
|
| 269 |
" -a missing_value,CMnight,o,b,255 ", |
|
| 270 |
" -a _FillValue,CMnight,o,b,255 ", |
|
| 271 |
" -a valid_range,CMnight,o,b,\"0,3\" ", |
|
| 272 |
" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ", |
|
| 273 |
ncfile,sep="")) |
|
| 274 |
#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep="")) |
|
| 360 | 275 |
|
| 361 |
|
|
| 362 |
### delete the temporary files |
|
| 363 |
unlink_.gislock() |
|
| 364 |
system(paste("rm -frR ",tf,sep=""))
|
|
| 365 |
|
|
| 366 | 276 |
|
| 367 | 277 |
## Confirm that the file has the correct attributes, otherwise delete it |
| 368 | 278 |
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
| ... | ... | |
| 373 | 283 |
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
| 374 | 284 |
file.remove(ncfile) |
| 375 | 285 |
} |
| 376 |
|
|
| 286 |
############ copy files to lou |
|
| 287 |
#if(platform=="pleiades"){
|
|
| 288 |
# archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
|
| 289 |
# system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
|
| 290 |
# system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
|
| 291 |
# file.remove(ncfile) |
|
| 292 |
# file.remove(paste(ncfile,".aux.xml",sep="")) |
|
| 293 |
#} |
|
| 294 |
|
|
| 295 |
|
|
| 296 |
### delete the temporary files |
|
| 297 |
# unlink_.gislock() |
|
| 298 |
# system(paste("rm -frR ",tempdir(),sep=""))
|
|
| 299 |
|
|
| 300 |
|
|
| 377 | 301 |
## print out some info |
| 378 |
print(paste(" ################################################################### Finished ",date,
|
|
| 379 |
"################################################################")) |
|
| 302 |
print(paste("####################### Finished ",tile,"-",date, "###################################"))
|
|
| 380 | 303 |
|
| 381 | 304 |
## delete old files |
| 382 |
system("cleartemp")
|
|
| 305 |
#system("cleartemp")
|
|
| 383 | 306 |
|
| 384 | 307 |
## quit |
| 385 | 308 |
q("no",status=0)
|
Also available in: Unified diff
updated to use new swtif