| 33 |
33 |
q(status=1);
|
| 34 |
34 |
}
|
| 35 |
35 |
|
|
36 |
testing=T
|
|
37 |
platform="pleiades"
|
| 36 |
38 |
|
| 37 |
39 |
## default date and tile to play with (will be overwritten below when running in batch)
|
| 38 |
|
date="20090129"
|
| 39 |
|
tile="h17v00"
|
| 40 |
|
platform="pleiades"
|
| 41 |
|
verbose=T
|
|
40 |
if(testing){
|
|
41 |
date="20090129"
|
|
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 |
|
|
|
47 |
if(!testing){
|
|
48 |
date=opta$date
|
|
49 |
tile=opta$tile
|
|
50 |
verbose=opta$verbose #print out extensive information for debugging?
|
|
51 |
}
|
| 48 |
52 |
## get year and doy from date
|
| 49 |
53 |
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
| 50 |
54 |
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
| ... | ... | |
| 55 |
59 |
## path to some executables
|
| 56 |
60 |
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
| 57 |
61 |
# swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif"
|
| 58 |
|
swtifpath="/nobackupp4/pvotava/software/heg/2.11/bin/swtif"
|
|
62 |
swtifpath="/nobackupp4/pvotava/software/heg/2.12/bin/swtif"
|
| 59 |
63 |
## path to swath database
|
| 60 |
64 |
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"
|
| 61 |
65 |
## specify working directory
|
| ... | ... | |
| 126 |
130 |
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk"))
|
| 127 |
131 |
|
| 128 |
132 |
|
| 129 |
|
############################################################################
|
| 130 |
|
############################################################################
|
| 131 |
|
### Use the HEG tool to grid all available swath data for this date-tile
|
| 132 |
|
for(file in swaths){
|
| 133 |
|
## Function to generate hegtool parameter file for multi-band HDF-EOS file
|
| 134 |
|
print(paste("Starting file",basename(file)))
|
| 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
|
| 143 |
|
## First write the parameter file (careful, heg is very finicky!)
|
| 144 |
|
hdr=paste("NUM_RUNS = ",length(km1vars$varid),"|MULTI_BAND_HDFEOS:",length(km1vars$varid),sep="")
|
| 145 |
|
grp=paste("
|
|
133 |
## define function that grids swaths
|
|
134 |
swtif<-function(file,var){
|
|
135 |
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path
|
|
136 |
## First write the parameter file (careful, heg is very finicky!)
|
|
137 |
hdr=paste("NUM_RUNS = 1")
|
|
138 |
grp=paste("
|
| 146 |
139 |
BEGIN
|
| 147 |
140 |
INPUT_FILENAME=",file,"
|
| 148 |
141 |
OBJECT_NAME=mod35
|
| 149 |
|
FIELD_NAME=",km1vars$variable,"|
|
| 150 |
|
BAND_NUMBER = 1
|
|
142 |
FIELD_NAME=",var$variable,"|
|
|
143 |
BAND_NUMBER = ",var$band,"
|
| 151 |
144 |
OUTPUT_PIXEL_SIZE_X=926.6
|
| 152 |
145 |
OUTPUT_PIXEL_SIZE_Y=926.6
|
| 153 |
146 |
# MODIS 1km Resolution
|
| 154 |
147 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
| 155 |
148 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
| 156 |
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",km1vars),"NN","CUBIC"),"
|
| 157 |
|
RESAMPLING_TYPE =NN
|
|
149 |
RESAMPLING_TYPE =",var$method,"
|
| 158 |
150 |
OUTPUT_PROJECTION_TYPE = SIN
|
| 159 |
151 |
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 )
|
| 160 |
152 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
| 161 |
153 |
ELLIPSOID_CODE = WGS84
|
| 162 |
154 |
OUTPUT_TYPE = HDFEOS
|
| 163 |
|
OUTPUT_FILENAME = ",km1outfile,"
|
|
155 |
OUTPUT_FILENAME = ",outfile,"
|
| 164 |
156 |
END
|
| 165 |
|
|
| 166 |
157 |
",sep="")
|
| 167 |
|
|
| 168 |
|
## if any remnants from previous runs remain, delete them
|
| 169 |
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
| 170 |
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
| 171 |
158 |
## write it to a file
|
| 172 |
159 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
| 173 |
160 |
## now run the swath2grid tool
|
| 174 |
161 |
## write the gridded file
|
| 175 |
162 |
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
|
|
163 |
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
|
164 |
}
|
| 205 |
165 |
|
| 206 |
|
",sep="")
|
|
166 |
## vars to grid
|
|
167 |
vars=as.data.frame(matrix(c(
|
|
168 |
"Cloud_Mask", "CM", "NN", 1,
|
|
169 |
"Quality_Assurance", "QA", "NN", 1,
|
|
170 |
"Solar_Zenith", "SolZen", "NN", 1,
|
|
171 |
"Sensor_Zenith", "SenZen", "CUBIC", 1
|
|
172 |
),
|
|
173 |
byrow=T,ncol=4,dimnames=list(1:4,c("variable","varid","method","band"))),stringsAsFactors=F)
|
| 207 |
174 |
|
| 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 |
175 |
|
| 214 |
|
print(paste("Finished gridding ", file," for tile ",tile))
|
|
176 |
############################################################################
|
|
177 |
############################################################################
|
|
178 |
### Use the HEG tool to grid all available swath data for this date-tile
|
|
179 |
for(file in swaths){
|
|
180 |
print(paste("Starting file",basename(file)))
|
|
181 |
## run swtif for each band
|
|
182 |
lapply(1:nrow(vars),function(i) swtif(file,vars[i,]))
|
| 215 |
183 |
} #end looping over swaths
|
| 216 |
184 |
|
| 217 |
|
########################
|
|
185 |
|
|
186 |
#############################################################################
|
|
187 |
## Check output
|
| 218 |
188 |
## confirm at least one file for this date is present. If not, quit.
|
| 219 |
189 |
outfiles=list.files(tempdir(),full=T,pattern=paste(basename(swaths),"$",sep="",collapse="|"))
|
| 220 |
190 |
if(!any(file.exists(outfiles))) {
|
| ... | ... | |
| 226 |
196 |
if(plot){
|
| 227 |
197 |
i=1
|
| 228 |
198 |
system(paste("gdalinfo ",swaths[1]))
|
| 229 |
|
d=brick(
|
| 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)
|
|
199 |
d=brick(lapply(outfiles,function(r) raster(r)))
|
|
200 |
plot(d)
|
| 236 |
201 |
}
|
| 237 |
202 |
#system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
| 238 |
203 |
|
| ... | ... | |
| 279 |
244 |
}
|
| 280 |
245 |
|
| 281 |
246 |
## Identify which files to process
|
| 282 |
|
tfs=unique(sub("km1_|km5_","",basename(outfiles)))
|
|
247 |
tfs=unique(sub("CM_|QA_|SenZen_|SolZen_","",basename(outfiles)))
|
| 283 |
248 |
#tfs=list.files(tempdir(),pattern="temp.*hdf")
|
| 284 |
249 |
nfs=length(tfs)
|
| 285 |
250 |
if(verbose) print(paste(nfs,"swaths available for processing"))
|
| 286 |
251 |
|
| 287 |
252 |
## loop through scenes and process QA flags
|
| 288 |
253 |
for(i in 1:nfs){
|
| 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(lat=melt(as.matrix(lat))[,3],lon=melt(as.matrix(lon))[,3],ID=1)
|
| 300 |
|
## crop to big bbox
|
| 301 |
|
coords=coords[coords$lat<tile_bb$lat_max+0.5&coords$lat>tile_bb$lat_min-0.5&
|
| 302 |
|
coords$lon>tile_bb$lon_min-0.5&coords$lon<tile_bb$lon_max+0.5,]
|
| 303 |
|
coordinates(coords)=c("lon","lat")
|
| 304 |
|
proj4string(coords)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
|
| 305 |
|
## project to sinusoidal
|
| 306 |
|
coords2=spTransform(coords,CRS(projection(td)))
|
| 307 |
|
writeOGR(coords2,dsn=".",layer=sub("[.]hdf","",basename(tswath)),driver="ESRI Shapefile",overwrite=T)
|
| 308 |
|
|
| 309 |
|
system(paste("gdal_grid -ot Byte -a count:radius1=10000:radius2=10000 ",
|
| 310 |
|
" -txe ",paste(bbox(td)[1,],sep="",collapse=" "),
|
| 311 |
|
" -tye ",paste(bbox(td)[2,],sep="",collapse=" "),
|
| 312 |
|
" -outsize ",paste(td@grid@cells.dim,collapse=" "),
|
| 313 |
|
" -l ",sub("[.]hdf","",basename(tswath))," ",
|
| 314 |
|
sub("[.]hdf",".shp",basename(tswath))," mask_",sub("[.]hdf",".tif",basename(tswath)),sep=""))
|
| 315 |
|
|
| 316 |
|
|
| 317 |
|
ps=rasterize(coords2,d2[[2]])
|
| 318 |
|
dist=distance(ps,edge=F)
|
| 319 |
|
dist2=dist<10000
|
| 320 |
|
|
| 321 |
|
## fit alpha hull to draw polygon around region with data
|
| 322 |
|
# ah=ahull(coordinates(coords2),alpha=100000)
|
| 323 |
|
# ah2=ah$x[ah$alpha.extremes,]
|
| 324 |
|
# ah2=ah$edges[,c("x1","y1")]
|
| 325 |
|
|
| 326 |
|
# pp = SpatialPolygons(list(Polygons(list(Polygon(coords[c(1:nrow(coords),1),])),1)))
|
| 327 |
|
# proj4string(pp)=projection(td)
|
| 328 |
|
|
| 329 |
|
|
| 330 |
|
plot(stack(lon,lat))
|
| 331 |
|
plot(coords,add=F);axis(1);axis(2)
|
| 332 |
|
plot(d2[[8]],add=F)
|
| 333 |
|
plot(coords2,add=T);axis(1);axis(2)
|
| 334 |
|
plot(ah,wpoints=F,add=T,col="red")
|
| 335 |
|
points(ah2,add=T)
|
| 336 |
|
|
| 337 |
|
glat=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLATITUDE=","",system(paste("gdalinfo ",tswath," | grep GRINGPOINTLATITUDE"),intern=T)),split=",")))
|
| 338 |
|
glon=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLONGITUDE=","",system(paste("gdalinfo ",tswath," | grep GRINGPOINTLONGITUDE"),intern=T)),split=",")))
|
| 339 |
|
bb=cbind(c(glon,glon[1]),c(glat,glat[1]))
|
| 340 |
|
pp = SpatialPolygons(list(Polygons(list(Polygon(bb)),1)))
|
| 341 |
|
proj4string(pp)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
|
| 342 |
|
pp2=spsample(as(pp,"SpatialLines"), 100, type="regular")
|
| 343 |
|
pp2=spTransform(pp2,CRS(projection(td)))
|
| 344 |
|
|
| 345 |
|
dt=projectRaster(d2[[1]],crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")
|
| 346 |
|
|
| 347 |
|
plot(pp,add=F,usePolypath = FALSE);axis(1);axis(2)#(sp.polygons(pp),usePolypath = FALSE)
|
| 348 |
|
plot(d2[[2]],add=F)
|
| 349 |
|
|
| 350 |
|
|
|
254 |
bfile=tfs[i]
|
|
255 |
## Read in the data from the HDFs
|
| 351 |
256 |
## Cloud Mask
|
| 352 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
|
|
257 |
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
| 353 |
258 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 354 |
|
## QA ## extract first bit to keep only "useful" values of cloud mask
|
| 355 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
|
|
259 |
## QA ## extract first bit to keep only "useful" values of cloud mask
|
|
260 |
execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),
|
| 356 |
261 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 357 |
|
## Sensor Zenith ## extract first bit to keep only "low angle" observations
|
| 358 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",km5file,"\":mod35:Sensor_Zenith",sep=""),
|
|
262 |
## Sensor Zenith ## extract first bit to keep only "low angle" observations
|
|
263 |
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
| 359 |
264 |
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 360 |
|
## Solar Zenith ## extract first bit to keep only "low angle" observations
|
| 361 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",km5file,"\":mod35:Solar_Zenith",sep=""),
|
|
265 |
|
|
266 |
## check for interpolation artefacts
|
|
267 |
# execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))
|
|
268 |
# execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))
|
|
269 |
# execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))
|
|
270 |
|
|
271 |
p=-50:50
|
|
272 |
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=""))
|
|
273 |
vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))
|
|
274 |
colnames(vals)=c("value","count")
|
|
275 |
vals$count=as.numeric(as.character(vals$count))
|
|
276 |
vals$value=as.numeric(as.character(vals$value))
|
|
277 |
vals=na.omit(vals)
|
|
278 |
vals$count[vals$value==1&vals$count>10]
|
|
279 |
#
|
|
280 |
#plot(p~value,data=vals)
|
|
281 |
# print(sum(vals$p[vals$p>.1]))
|
|
282 |
|
|
283 |
## TODO: use this to flag problem swaths?
|
|
284 |
## Solar Zenith ## extract first bit to keep only "low angle" observations
|
|
285 |
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
| 362 |
286 |
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 363 |
|
|
|
287 |
## produce the summaries
|
| 364 |
288 |
system(paste("r.mapcalc <<EOF
|
| 365 |
289 |
CM_fill_",i," = if(isnull(CM1_",i,"),1,0)
|
| 366 |
290 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0)
|
| 367 |
|
SZ_low_",i," = if(SZ_",i,"<6000,1,0)
|
|
291 |
SZ_low_",i," = if(SZ_",i,"_clump==0&SZ_",i,"<6000,1,0)
|
| 368 |
292 |
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0)
|
| 369 |
293 |
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
| 370 |
|
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999)
|
| 371 |
|
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",-9999)
|
| 372 |
|
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",-9999)
|
|
294 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null())
|
|
295 |
SZday_",i," = if(SZ_",i,"_clump==0&CM_dayflag_",i,"==1,SZ_",i,",null())
|
|
296 |
SZnight_",i," = if(SZ_",i,"_clump==0&CM_dayflag_",i,"==0,SZ_",i,",null())
|
|
297 |
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())
|
|
298 |
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())
|
| 373 |
299 |
EOF",sep=""))
|
| 374 |
300 |
|
| 375 |
301 |
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
| 376 |
|
# 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))))
|
| 377 |
|
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))
|
| 378 |
|
|
| 379 |
|
## set null values
|
| 380 |
|
execGRASS("r.null",map=paste("CMday_",i,sep=""),setnull="-9999")
|
| 381 |
|
execGRASS("r.null",map=paste("CMnight_",i,sep=""),setnull="-9999")
|
|
302 |
# 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))))
|
|
303 |
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))
|
| 382 |
304 |
|
| 383 |
305 |
drawplot=F
|
| 384 |
306 |
if(drawplot){
|
| 385 |
307 |
d2=stack(
|
| 386 |
|
raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
| 387 |
|
raster(readRAST6(paste("CM1_",i,sep=""))),
|
| 388 |
|
raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
| 389 |
|
raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
308 |
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
|
309 |
# raster(readRAST6(paste("CM1_",i,sep=""))),
|
|
310 |
# raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
311 |
# raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
312 |
# raster(readRAST6(paste("CMday_",i,sep=""))),
|
| 390 |
313 |
raster(readRAST6(paste("CMnight_",i,sep=""))),
|
| 391 |
|
raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
| 392 |
|
raster(readRAST6(paste("SoZ_",i,sep=""))),
|
| 393 |
|
raster(readRAST6(paste("SZ_",i,sep="")))
|
|
314 |
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
|
315 |
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
|
316 |
raster(readRAST6(paste("SZ_",i,sep=""))),
|
|
317 |
raster(readRAST6(paste("SZ_",i,"_clump",sep="")))
|
| 394 |
318 |
)
|
| 395 |
|
plot(d2[[2]],add=F)
|
| 396 |
|
points(coords2,pch=16,cex=.2,add=T)
|
| 397 |
|
plot(pp,add=F,usePolypath = FALSE)#(sp.polygons(pp),usePolypath = FALSE)
|
| 398 |
|
|
|
319 |
plot(d2,add=F)
|
| 399 |
320 |
}
|
| 400 |
321 |
|
| 401 |
322 |
|
| ... | ... | |
| 403 |
324 |
|
| 404 |
325 |
## select lowest view angle
|
| 405 |
326 |
## use r.series to find minimum
|
| 406 |
|
system("r.series input=`g.mlist rast pat=\"SZ_[0-9]*$\" sep=,` output=SZ_min method=min_raster")
|
|
327 |
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
|
328 |
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
| 407 |
329 |
|
| 408 |
330 |
## select cloud observation with lowest sensor zenith for day and night
|
| 409 |
331 |
system(
|
| 410 |
332 |
paste("r.mapcalc <<EOF
|
| 411 |
|
CMday_daily=",paste(paste("if((SZ_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
| 412 |
|
CMnight_daily=",paste(paste("if((SZ_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
|
333 |
CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
|
334 |
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
| 413 |
335 |
))
|
| 414 |
336 |
|
| 415 |
|
d=brick(lapply(1:nfs,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
| 416 |
|
levelplot(d)
|
| 417 |
|
|
| 418 |
|
plot(raster(readRAST6("CMnight_daily")))
|
| 419 |
|
|
| 420 |
|
#### Now generate daily minimum p(clear)
|
| 421 |
|
system(paste("r.mapcalc <<EOF
|
| 422 |
|
CMday_daily=int((min(",paste("if(isnull(CMday_",1:nfs,"),9999,CMday_",1:nfs,")",sep="",collapse=","),")))
|
| 423 |
|
CMnight_daily=int((min(",paste("if(isnull(CMnight_",1:nfs,"),9999,CMnight_",1:nfs,")",sep="",collapse=","),")))
|
| 424 |
|
EOF",sep=""))
|
| 425 |
|
|
|
337 |
if(plot){
|
|
338 |
sz1=brick(lapply(1:nfs,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))
|
|
339 |
sz_clump=brick(lapply(1:nfs,function(i) raster(readRAST6(paste("SZ_",i,"_clump",sep="")))))
|
|
340 |
d=brick(lapply(1:nfs,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
341 |
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
|
342 |
library(rasterVis)
|
|
343 |
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))
|
|
344 |
levelplot(sz_clump)
|
|
345 |
levelplot(d)
|
|
346 |
levelplot(d2)
|
|
347 |
}
|
| 426 |
348 |
|
| 427 |
|
## reset null values
|
| 428 |
|
execGRASS("r.null",map="CMday_daily",setnull="9999")
|
| 429 |
|
execGRASS("r.null",map="CMnight_daily",setnull="9999")
|
| 430 |
349 |
|
| 431 |
350 |
### Write the files to a netcdf file
|
| 432 |
351 |
## create image group to facilitate export as multiband netcdf
|
| ... | ... | |
| 496 |
415 |
|
| 497 |
416 |
|
| 498 |
417 |
## print out some info
|
| 499 |
|
print(paste(" ################################################################### Finished ",date,
|
|
418 |
print(paste(" ################################################################### Finished ",tile,"-",date,
|
| 500 |
419 |
"################################################################"))
|
| 501 |
420 |
|
| 502 |
421 |
## delete old files
|
added function to flag bad pixels introduced by bug in HEG software