Quality Assessment of the Processing of “Vol 3 Config 2” (TTC3322)
Requirements and Settings
root <- "/media/alobo"
#root <- "/Volumes"
options(digits.secs=3)
require(rgdal)
require(raster)
source(paste(root,"/LACIE500/Rutils/RLRS/ima6b.R",sep=""))
rdir <- paste(root,"/LACIE500/LIFEN/castello20140708/vol3_config2",sep="")
w <- 80
imafname <- "TTC3522.tif"
indfname <- "TTC3522ConIndices.tif"1. Read images
# Mapa de la configuración ( de máster a banda 5):
# R780
# R530
# R570
# R670
# R710
# R730 1.1 Convert PW2 tif to multiband geotif
ima6b(imasdir=rdir,nomin=imafname,imasdirout=rdir, nomout="TTC3522v2.tif",
b16=TRUE,mos=F,orden=c(2:6,1),
flip=FALSE, naflag=0, plotout=TRUE, namebands=NA,geo=TRUE)
imav3c2 <- brick(file.path(rdir,"TTC3522v2.tif"))
imav3c2[imav3c2==0] <- NA
names(imav3c2) <- c("R530", "R570","R670","R710","R730","R780")
options(width = 180)
summary(imav3c2)## TTC3522v2.1 TTC3522v2.2 TTC3522v2.3 TTC3522v2.4 TTC3522v2.5 TTC3522v2.6
## Min. 4 4 4 4 4 113
## 1st Qu. 241 361 441 349 356 414
## Median 306 489 596 456 429 494
## 3rd Qu. 394 626 761 560 522 595
## Max. 633 1014 1022 821 763 1020
## NA's 12297 9579 12188 13109 20961 01.2. Reflectance image
“Bands” 2 to 7 in “TTC3522ConIndices.tif” (Hemav format)
indv3c2 <- brick(file.path(rdir,indfname))
imareflv3c2 <- indv3c2[[c(3:7,2)]]
names(imareflv3c2) <- c("R530", "R570","R670","R710","R730","R780")
imareflv3c2[imareflv3c2==0] <- NA
summary(imareflv3c2)## R530 R570 R670 R710 R730 R780
## Min. 4 4 53 4 77 113
## 1st Qu. 241 361 442 349 356 415
## Median 306 489 596 457 429 495
## 3rd Qu. 394 628 761 561 522 596
## Max. 633 1014 1022 821 763 1020
## NA's 12297 18485 20885 23958 29687 296872 Compare uncalibrated and calibrated images
2.1 CIR images
With axes=TRUE
par(mfrow=c(1,2))
plotRGB(imav3c2,r=6,g=3,b=2,stretch="lin",main="Uncalbrated CIR",axes=TRUE)
plotRGB(imareflv3c2,r=6,g=3,b=2,stretch="lin",main="Calibrated CIR",axes=TRUE)
With axes=FALSE
par(mfrow=c(1,2))
plotRGB(imav3c2,r=6,g=3,b=2,stretch="lin",main="Uncalbrated CIR",axes=FALSE)
plotRGB(imareflv3c2,r=6,g=3,b=2,stretch="lin",main="Calibrated CIR",axes=FALSE)
Uncalibrated and calibrated images are identical!!* ### 2.2 Inspect the same random set of points
rndp1 <- sampleRandom(imav3c2,50,rowcol=TRUE)
rndp2 <- cbind(rndp1[,1:2],imareflv3c2[rndp1[,1:2]])
#rndp <- stack(imav3c2,imareflv3c2)[rndp1[,1:2]]
rndpv3c2 <- rbind(rndp1,rndp2)
rndpv3c2 <- cbind(1,rndpv3c2)
rndpv3c2[51:100,1] <- 2
rm(rndp1,rndp2)
#dev.off()
par(mfrow=c(1,1))
matplot(x=c(530,570,670,710,730,780),t(rndpv3c2[,-(1:3)]),type="b",pch=18,
xlab="Wavelength (nm)",ylab="Reflectance",col=rndpv3c2[,1])
plot(x=c(530,570,670,710,730,780), y=rndpv3c2[1,-(1:3)],type="b")
points(x=c(530,570,670,710,730,780), y=rndpv3c2[51,-(1:3)],type="b",col="red")
Spectra are identical in Calibrated and not calibrated???
#dev.off()
par(mfrow=c(2,1))
plot(imav3c2[[1]],imareflv3c2[[1]],grid=TRUE,xlab="Not Calibrated",ylab="Calibrated",main="R450",asp=1)
plot(imav3c2[[6]],imareflv3c2[[6]],grid=TRUE,xlab="Not Calibrated",ylab="Calibrated",main="R800",asp=1)
Images have not been calibrated ==> Indices will be wrong ## 3. Indices ### 3.1 Hemav indices
# Indices calculados:
# Banda número 1: Banda con Stretching
# Banda número 8: NDVI
# Banda número 9: DCNI
# Banda número 10: PRI
indv3c2 <- brick(file.path(rdir,indfname))Select indices only
indv3c2 <- indv3c2[[8:10]]
names(indv3c2) <- c("ndvi","dcni","pri")3.2 Calculate indices from the (supposedly) calibrated image
NDVI Definition: (R800 − R670) / (R800 + R670)
Note we use 780 instead of 800
imareflv3c2ndvi <- (imareflv3c2[["R780"]] - imareflv3c2[["R670"]])/(imareflv3c2[["R780"]] + imareflv3c2[["R670"]])
summary(imareflv3c2ndvi)## layer
## Min. -4.740e-01
## 1st Qu. -1.517e-01
## Median -1.230e-01
## 3rd Qu. -6.888e-02
## Max. 8.528e-01
## NA's 2.969e+04NDVI values are too high for this scene
delme <- writeRaster(imareflv3c2ndvi,file="imareflv3c2ndvi780",overwrite=TRUE, format="GTiff",NAflag=-9999)DCNI Definition: ((R720 − R700) / (R700 − R670)) / (R720 − R670 + 0.03)
Note we use 730 and 710 instead of 720 and 700
imareflv3c2dcni <- ((imareflv3c2[["R730"]]-imareflv3c2[["R710"]])/(imareflv3c2[["R710"]]-imareflv3c2[["R670"]]))/
(imareflv3c2[["R730"]]-imareflv3c2[["R670"]] +0.03)
summary(imareflv3c2dcni)## layer
## Min. -3.333e+01
## 1st Qu. -1.693e-03
## Median -8.133e-04
## 3rd Qu. -1.057e-04
## Max. Inf
## NA's 2.970e+04delme <- writeRaster(imareflv3c2dcni,file="imareflv3c2dcni2",overwrite=TRUE, format="GTiff",NAflag=-9999)Note: weird values, perhaps because uncalibrated image
PRI570 Definition: (R570 − R531) / (R570 + R531)
Note we use 530 instead of 531
imareflv3c2pri <- (imareflv3c2[["R570"]] - imareflv3c2[["R530"]])/(imareflv3c2[["R570"]] + imareflv3c2[["R530"]])
summary(imareflv3c2pri)## layer
## Min. -0.9821
## 1st Qu. 0.1948
## Median 0.2213
## 3rd Qu. 0.2390
## Max. 0.9903
## NA's 18485.0000delme <- writeRaster(imareflv3c2pri,file="imareflv3c2pri",overwrite=TRUE, format="GTiff",NAflag=-9999)3.3 Compare NDVI
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2ndvi, legend=FALSE, axes=FALSE, main="ALOBO NDVI")
plot(indv3c2[["ndvi"]], legend=FALSE, axes=FALSE, main="Hemav NDVI")
par(mfrow=c(1,2))
hist(imareflv3c2ndvi,main="ALOBO",xlab="NDVI")
hist(indv3c2[["ndvi"]], main="Hemav",xlab="NDVI*1000")
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2ndvi,indv3c2[["ndvi"]],xlab="NDVI",ylab="Hemav_NDVI",pch=20)
#dev.off()Note they are different!!! ### 3.4 Compare DCNI
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2dcni, legend=FALSE, axes=FALSE, main="ALOBO DCNI")
plot(indv3c2[["dcni"]], legend=FALSE, axes=FALSE, main="Hemav DCNI*1000")
par(mfrow=c(1,2))
hist(imareflv3c2dcni,main="ALOBO",xlab="ALOBO dcni")
hist(indv3c2[["dcni"]], main="Hemav",xlab="Hemav dcni*1000")
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2dcni,indv3c2[["dcni"]],xlab="ALOBO dcni",ylab="Hemav_dcni*1000",pch=20)
#dev.off()Weird result (ALOBO), probably because uncalibrated images ### 3.5 Compare PRI
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2pri, legend=FALSE, axes=FALSE, main="ALOBO PRI")
plot(indv3c2[["pri"]], legend=FALSE, axes=FALSE, main="Hemav PRI")
par(mfrow=c(1,2))
hist(imareflv3c2pri,main="ALOBO",xlab="PRI")
hist(indv3c2[["pri"]], main="Hemav",xlab="PRI*1000")
#dev.off()
par(mfrow=c(1,1))
plot(imareflv3c2pri,indv3c2[["pri"]],xlab="PRI",ylab="Hemav_PRI",pch=20)
Conclusions
- Images in “bands” 2 to 7 in TTC3322ConIndices.tif are not calibrated
- Indices must be computed with calibrated images
- DCNI values are weird. Probably because of using uncalibrated images
- Al 3 indices in “bands” 8 to 10 in TTC3322ConIndices.tif are very different than my calculations
- I do not have the saturating NDVI problem