Revision a981dca4
Added by Adam Wilson over 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="20090129"
|
|
| 39 |
#tile="h17v00"
|
|
| 38 |
date="20090129" |
|
| 39 |
tile="h17v00" |
|
| 40 | 40 |
platform="pleiades" |
| 41 | 41 |
verbose=T |
| 42 | 42 |
|
| ... | ... | |
| 44 | 44 |
date=opta$date |
| 45 | 45 |
tile=opta$tile |
| 46 | 46 |
verbose=opta$verbose #print out extensive information for debugging? |
| 47 |
|
|
| 47 | 48 |
## get year and doy from date |
| 48 | 49 |
year=format(as.Date(date,"%Y%m%d"),"%Y") |
| 49 | 50 |
doy=format(as.Date(date,"%Y%m%d"),"%j") |
| ... | ... | |
| 53 | 54 |
datadir=paste("/nobackupp4/datapool/modis/MOD35_L2.006/",year,"/",doy,"/",sep="")
|
| 54 | 55 |
## path to some executables |
| 55 | 56 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/" |
| 56 |
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif" |
|
| 57 |
swtifpath="/nobackupp4/pvotava/software/heg/bin/swtif" |
|
| 58 |
## path to swath database |
|
| 59 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
|
| 60 |
## specify working directory |
|
| 61 |
outdir=paste("/nobackupp1/awilso10/mod35/daily/",tile,"/",sep="") #directory for separate daily files
|
|
| 62 |
basedir="/nobackupp1/awilso10/mod35/" #directory to hold files temporarily before transferring to lou |
|
| 63 |
setwd(tempdir()) |
|
| 64 |
## grass database |
|
| 65 |
gisBase="/u/armichae/pr/grass-6.4.2/" |
|
| 66 |
## path to MOD11A1 file for this tile to align grid/extent |
|
| 67 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A2.005/2009.01.01",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
|
| 68 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_8Day_1km_LST:LST_Day_1km",sep=""))
|
|
| 69 |
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 " |
|
| 70 |
} |
|
| 57 |
# swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif" |
|
| 58 |
swtifpath="/nobackupp4/pvotava/software/heg/2.11/bin/swtif" |
|
| 71 | 59 |
## path to swath database |
| 72 | 60 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db" |
| 73 | 61 |
## specify working directory |
| ... | ... | |
| 114 | 102 |
upleft=paste(min(90,tile_bb$lat_max+0.5),max(-180,tile_bb$lon_min-0.5)) #northwest corner |
| 115 | 103 |
lowright=paste(max(-90,tile_bb$lat_min-0.5),min(180,tile_bb$lon_max+0.5)) #southeast corner |
| 116 | 104 |
|
| 117 |
## vars to process |
|
| 118 |
vars=as.data.frame(matrix(c( |
|
| 119 |
"Cloud_Mask", "CM", |
|
| 120 |
"Quality_Assurance", "QA"), |
|
| 121 |
# "Solar_Azimuth", "SA", |
|
| 122 |
# "Sensor_Zenith", "SZ"), |
|
| 123 |
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
|
| 124 |
|
|
| 125 | 105 |
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end. |
| 126 | 106 |
finalvars=c("CMday","CMnight")
|
| 127 | 107 |
|
| ... | ... | |
| 152 | 132 |
for(file in swaths){
|
| 153 | 133 |
## Function to generate hegtool parameter file for multi-band HDF-EOS file |
| 154 | 134 |
print(paste("Starting file",basename(file)))
|
| 155 |
outfile=paste(tempdir(),"/",basename(file),sep="") |
|
| 156 |
### Get 1km data |
|
| 135 |
|
|
| 136 |
## vars to process |
|
| 137 |
km1vars=as.data.frame(matrix(c( |
|
| 138 |
"Cloud_Mask", "CM", |
|
| 139 |
"Quality_Assurance", "QA"), |
|
| 140 |
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
|
| 141 |
km1outfile=paste(tempdir(),"/km1_",basename(file),sep="") |
|
| 142 |
## Get 1km data |
|
| 157 | 143 |
## First write the parameter file (careful, heg is very finicky!) |
| 158 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 144 |
hdr=paste("NUM_RUNS = ",length(km1vars$varid),"|MULTI_BAND_HDFEOS:",length(km1vars$varid),sep="")
|
|
| 159 | 145 |
grp=paste("
|
| 160 | 146 |
BEGIN |
| 161 | 147 |
INPUT_FILENAME=",file," |
| 162 | 148 |
OBJECT_NAME=mod35 |
| 163 |
FIELD_NAME=",vars$variable,"| |
|
| 149 |
FIELD_NAME=",km1vars$variable,"|
|
|
| 164 | 150 |
BAND_NUMBER = 1 |
| 165 | 151 |
OUTPUT_PIXEL_SIZE_X=926.6 |
| 166 | 152 |
OUTPUT_PIXEL_SIZE_Y=926.6 |
| 167 | 153 |
# MODIS 1km Resolution |
| 168 | 154 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," ) |
| 169 | 155 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," ) |
| 170 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
|
| 156 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",km1vars),"NN","CUBIC"),"
|
|
| 171 | 157 |
RESAMPLING_TYPE =NN |
| 172 | 158 |
OUTPUT_PROJECTION_TYPE = SIN |
| 173 | 159 |
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 ) |
| 174 | 160 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
| 175 | 161 |
ELLIPSOID_CODE = WGS84 |
| 176 | 162 |
OUTPUT_TYPE = HDFEOS |
| 177 |
OUTPUT_FILENAME = ",outfile," |
|
| 163 |
OUTPUT_FILENAME = ",km1outfile,"
|
|
| 178 | 164 |
END |
| 179 | 165 |
|
| 180 | 166 |
",sep="") |
| ... | ... | |
| 187 | 173 |
## now run the swath2grid tool |
| 188 | 174 |
## write the gridded file |
| 189 | 175 |
system(paste(swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -tmpLatLondir ",tempdir(),sep=""),intern=F,ignore.stderr=F) |
| 176 |
################# |
|
| 177 |
## 5km quality data |
|
| 178 |
km5vars=as.data.frame(matrix(c( |
|
| 179 |
"Solar_Zenith", "SolZen", |
|
| 180 |
"Sensor_Zenith", "SenZen"), |
|
| 181 |
byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
|
|
| 182 |
km5outfile=paste(tempdir(),"/km5_",basename(file),sep="") |
|
| 183 |
|
|
| 184 |
hdr=paste("NUM_RUNS = ",length(km5vars$varid),"|MULTI_BAND_HDFEOS:",length(km5vars$varid),sep="")
|
|
| 185 |
grp=paste("
|
|
| 186 |
BEGIN |
|
| 187 |
INPUT_FILENAME=",file," |
|
| 188 |
OBJECT_NAME=mod35 |
|
| 189 |
FIELD_NAME=",km5vars$variable,"| |
|
| 190 |
BAND_NUMBER = 1 |
|
| 191 |
OUTPUT_PIXEL_SIZE_X=926.6 |
|
| 192 |
OUTPUT_PIXEL_SIZE_Y=926.6 |
|
| 193 |
# MODIS 1km Resolution |
|
| 194 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," ) |
|
| 195 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," ) |
|
| 196 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",km1vars),"NN","CUBIC"),"
|
|
| 197 |
RESAMPLING_TYPE =CUBIC |
|
| 198 |
OUTPUT_PROJECTION_TYPE = SIN |
|
| 199 |
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 ) |
|
| 200 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
|
| 201 |
ELLIPSOID_CODE = WGS84 |
|
| 202 |
OUTPUT_TYPE = HDFEOS |
|
| 203 |
OUTPUT_FILENAME = ",km5outfile," |
|
| 204 |
END |
|
| 205 |
|
|
| 206 |
",sep="") |
|
| 207 |
|
|
| 208 |
## write it to a file |
|
| 209 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
|
| 210 |
## now run the swath2grid tool |
|
| 211 |
## write the gridded file |
|
| 212 |
system(paste(swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -tmpLatLondir ",tempdir(),sep=""),intern=F,ignore.stderr=F) |
|
| 213 |
|
|
| 190 | 214 |
print(paste("Finished gridding ", file," for tile ",tile))
|
| 191 | 215 |
} #end looping over swaths |
| 192 | 216 |
|
| 193 | 217 |
######################## |
| 194 | 218 |
## confirm at least one file for this date is present. If not, quit. |
| 195 |
outfiles=paste(tempdir(),"/",basename(swaths),sep="")
|
|
| 219 |
outfiles=list.files(tempdir(),full=T,pattern=paste(basename(swaths),"$",sep="",collapse="|"))
|
|
| 196 | 220 |
if(!any(file.exists(outfiles))) {
|
| 197 | 221 |
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
| 198 | 222 |
q("no",status=0)
|
| ... | ... | |
| 201 | 225 |
plot=F |
| 202 | 226 |
if(plot){
|
| 203 | 227 |
i=1 |
| 204 |
GDALinfo(outfiles[i])
|
|
| 228 |
system(paste("gdalinfo ",swaths[1]))
|
|
| 205 | 229 |
d=brick( |
| 206 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[i],"\":mod35:Cloud_Mask_0",sep="")),
|
|
| 207 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[i],"\":mod35:Quality_Assurance_0",sep="")),
|
|
| 208 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[i],"\":mod35:Sensor_Zenith",sep="")),
|
|
| 209 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[i],"\":mod35:Solar_Azimuth",sep="")))
|
|
| 210 |
plot(d) |
|
| 230 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[1],"\":mod35:Cloud_Mask_0",sep="")),
|
|
| 231 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[1],"\":mod35:Quality_Assurance_0",sep="")),
|
|
| 232 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[26],"\":mod35:Sensor_Zenith",sep="")),
|
|
| 233 |
raster(paste("HDF4_EOS:EOS_GRID:\"",outfiles[26],"\":mod35:Solar_Zenith",sep=""))
|
|
| 234 |
) |
|
| 235 |
plot(d) |
|
| 211 | 236 |
} |
| 212 |
# TODO: import sensor zenith separately and mask out bad pixels |
|
| 213 |
system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
|
| 237 |
#system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
|
| 214 | 238 |
|
| 215 | 239 |
##################################################### |
| 216 | 240 |
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths |
| ... | ... | |
| 237 | 261 |
|
| 238 | 262 |
## set up temporary grass instance for this PID |
| 239 | 263 |
if(verbose) print(paste("Set up temporary grass session in",tf))
|
| 240 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 264 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod35",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
|
| 241 | 265 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 242 | 266 |
|
| 243 | 267 |
## Define region by importing one MOD11A1 raster. |
| ... | ... | |
| 253 | 277 |
output="modisgrid",flags=c("overwrite","o"))
|
| 254 | 278 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
| 255 | 279 |
} |
| 280 |
|
|
| 256 | 281 |
## Identify which files to process |
| 257 |
tfs=basename(swaths) |
|
| 258 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
|
| 259 |
tfs=tfs[tfs%in%list.files(tempdir())] |
|
| 282 |
tfs=unique(sub("km1_|km5_","",basename(outfiles)))
|
|
| 260 | 283 |
#tfs=list.files(tempdir(),pattern="temp.*hdf") |
| 261 | 284 |
nfs=length(tfs) |
| 262 | 285 |
if(verbose) print(paste(nfs,"swaths available for processing")) |
| 263 | 286 |
|
| 264 | 287 |
## loop through scenes and process QA flags |
| 265 | 288 |
for(i in 1:nfs){
|
| 266 |
file=paste(tempdir(),"/",tfs[i],sep="") |
|
| 289 |
file=paste(tempdir(),"/km1_",tfs[i],sep="") |
|
| 290 |
km5file=paste(tempdir(),"/km5_",tfs[i],sep="") |
|
| 291 |
|
|
| 292 |
## get GRING coordinates of swath to crop raster |
|
| 293 |
tswath=swaths[grep(tfs[i],swaths)] |
|
| 294 |
system(paste("gdalinfo ",tswath))
|
|
| 295 |
lat=raster(paste("HDF4_SDS:UNKNOWN:\"",tswath,"\":0",sep=""))
|
|
| 296 |
lon=raster(paste("HDF4_SDS:UNKNOWN:\"",tswath,"\":1",sep=""))
|
|
| 297 |
## HEG Tool reprojection results in large areas of sprious values (regions outside data areas are filled in using the interpolation method) |
|
| 298 |
## need to crop the resulting projected data to eliminate these areas |
|
| 299 |
coords=cbind.data.frame(melt(as.matrix(lat))[,1:3],lon=melt(as.matrix(lon))[,3]) |
|
| 300 |
coords=coords[coords$X1%in%range(coords$X1)|coords$X2%in%range(coords$X2),4:3];colnames(coords)=c("lon","lat")
|
|
| 301 |
#coords=cbind.data.frame(lat=melt(as.matrix(lat))[,3],lon=melt(as.matrix(lon))[,3]) |
|
| 302 |
## crop to big bbox |
|
| 303 |
# coords=coords[coords$lat<tile_bb$lat_max+0.5&coords$lat>tile_bb$lat_min-0.5& |
|
| 304 |
# coords$lon>tile_bb$lon_min-0.5&coords$lon<tile_bb$lon_max+0.5,] |
|
| 305 |
# coordinates(coords)=c("lon","lat")
|
|
| 306 |
# proj4string(coords)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs " |
|
| 307 |
## project to sinusoidal |
|
| 308 |
coords2=spTransform(coords,CRS(projection(td))) |
|
| 309 |
## fit alpha hull to draw polygon around region with data |
|
| 310 |
ah=ahull(coordinates(coords2),alpha=100000) |
|
| 311 |
ah2=ah$x[ah$alpha.extremes,] |
|
| 312 |
ah2=ah$edges[,c("x1","y1")]
|
|
| 313 |
|
|
| 314 |
pp = SpatialPolygons(list(Polygons(list(Polygon(coords[c(1:nrow(coords),1),])),1))) |
|
| 315 |
proj4string(pp)=projection(td) |
|
| 316 |
|
|
| 317 |
plot(coords,add=F);axis(1);axis(2) |
|
| 318 |
plot(d2[[2]],add=F) |
|
| 319 |
plot(coords2,add=T);axis(1);axis(2) |
|
| 320 |
plot(ah,wpoints=F,add=T,col="red") |
|
| 321 |
points(ah2,add=T) |
|
| 322 |
|
|
| 323 |
glat=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLATITUDE=","",system(paste("gdalinfo ",tswath," | grep GRINGPOINTLATITUDE"),intern=T)),split=",")))
|
|
| 324 |
glon=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLONGITUDE=","",system(paste("gdalinfo ",tswath," | grep GRINGPOINTLONGITUDE"),intern=T)),split=",")))
|
|
| 325 |
bb=cbind(c(glon,glon[1]),c(glat,glat[1])) |
|
| 326 |
pp = SpatialPolygons(list(Polygons(list(Polygon(bb)),1))) |
|
| 327 |
proj4string(pp)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs " |
|
| 328 |
pp2=spsample(as(pp,"SpatialLines"), 100, type="regular") |
|
| 329 |
pp2=spTransform(pp2,CRS(projection(td))) |
|
| 330 |
|
|
| 331 |
dt=projectRaster(d2[[1]],crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ") |
|
| 332 |
|
|
| 333 |
plot(pp,add=F,usePolypath = FALSE);axis(1);axis(2)#(sp.polygons(pp),usePolypath = FALSE) |
|
| 334 |
plot(d2[[2]],add=F) |
|
| 335 |
|
|
| 336 |
|
|
| 267 | 337 |
## Cloud Mask |
| 268 | 338 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
|
| 269 | 339 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 270 | 340 |
## QA ## extract first bit to keep only "useful" values of cloud mask |
| 271 | 341 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
|
| 272 | 342 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 273 |
## extract first two bits of cloud mask |
|
| 274 |
system(paste("r.mapcalc <<EOF
|
|
| 343 |
## Sensor Zenith ## extract first bit to keep only "low angle" observations |
|
| 344 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",km5file,"\":mod35:Sensor_Zenith",sep=""),
|
|
| 345 |
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 346 |
## Solar Zenith ## extract first bit to keep only "low angle" observations |
|
| 347 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",km5file,"\":mod35:Solar_Zenith",sep=""),
|
|
| 348 |
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 349 |
|
|
| 350 |
system(paste("r.mapcalc <<EOF
|
|
| 351 |
CM_fill_",i," = if(isnull(CM1_",i,"),1,0) |
|
| 275 | 352 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0) |
| 276 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999) |
|
| 353 |
SZ_low_",i," = if(SZ_",i,"<6000,1,0) |
|
| 354 |
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0) |
|
| 277 | 355 |
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0) |
| 278 |
CMday_",i," = if(QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",-9999) |
|
| 279 |
CMnight_",i," = if(QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",-9999) |
|
| 356 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999) |
|
| 357 |
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",-9999) |
|
| 358 |
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",-9999) |
|
| 280 | 359 |
EOF",sep="")) |
| 281 |
#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)))) |
|
| 282 |
#CM_path_",i," = ((CM1_",i," / 2^6) % 2^2) |
|
| 283 | 360 |
|
| 284 |
# execGRASS("r.null",map=paste("Pclearday_",i,sep=""),setnull="-9999")
|
|
| 285 |
# execGRASS("r.null",map=paste("Pclearnight_",i,sep=""),setnull="-9999")
|
|
| 361 |
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0) |
|
| 362 |
# 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)))) |
|
| 363 |
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4))) |
|
| 286 | 364 |
|
| 365 |
## set null values |
|
| 287 | 366 |
execGRASS("r.null",map=paste("CMday_",i,sep=""),setnull="-9999")
|
| 288 | 367 |
execGRASS("r.null",map=paste("CMnight_",i,sep=""),setnull="-9999")
|
| 289 |
|
|
| 290 |
|
|
| 368 |
|
|
| 369 |
drawplot=F |
|
| 370 |
if(drawplot){
|
|
| 371 |
d2=stack( |
|
| 372 |
raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
|
| 373 |
raster(readRAST6(paste("CM1_",i,sep=""))),
|
|
| 374 |
raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
| 375 |
raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
| 376 |
raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
| 377 |
raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
|
| 378 |
raster(readRAST6(paste("SoZ_",i,sep=""))),
|
|
| 379 |
raster(readRAST6(paste("SZ_",i,sep="")))
|
|
| 380 |
) |
|
| 381 |
plot(d2[[2]],add=F) |
|
| 382 |
plot(coords2,pch=16,cex=.2,add=T) |
|
| 383 |
plot(pp,add=F,usePolypath = FALSE)#(sp.polygons(pp),usePolypath = FALSE) |
|
| 384 |
} |
|
| 385 |
|
|
| 386 |
|
|
| 291 | 387 |
} #end loop through sub daily files |
| 292 | 388 |
|
| 389 |
## select lowest view angle |
|
| 390 |
## use r.series to find minimum |
|
| 391 |
system("r.series input=`g.mlist rast pat=\"SZ_[0-9]*$\" sep=,` output=SZ_min method=min_raster")
|
|
| 392 |
|
|
| 393 |
## select cloud observation with lowest sensor zenith for day and night |
|
| 394 |
system( |
|
| 395 |
paste("r.mapcalc <<EOF
|
|
| 396 |
CMday_daily=",paste(paste("if((SZ_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
|
| 397 |
CMnight_daily=",paste(paste("if((SZ_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
|
| 398 |
)) |
|
| 399 |
|
|
| 400 |
d=brick(lapply(1:nfs,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
| 401 |
levelplot(d) |
|
| 402 |
|
|
| 403 |
plot(raster(readRAST6("CMnight_daily")))
|
|
| 404 |
|
|
| 293 | 405 |
#### Now generate daily minimum p(clear) |
| 294 | 406 |
system(paste("r.mapcalc <<EOF
|
| 295 | 407 |
CMday_daily=int((min(",paste("if(isnull(CMday_",1:nfs,"),9999,CMday_",1:nfs,")",sep="",collapse=","),")))
|
Also available in: Unified diff
Adding sensor zenith to quality calculation and improving the filter used to select daily data