################################################### ### chunk number 1: setuphide ################################################### options(width=75) # load the landsat package library("landsat") # import July images data("july1") data("july2") data("july3") data("july4") data("july5") data("july61") data("july62") data("july7") # import November images data("nov1") data("nov2") data("nov3") data("nov4") data("nov5") data("nov61") data("nov62") data("nov7") # import digital elevation model data("dem") # add a modified image() function imageSGDF <- function(x, col = heat.colors(12), zlim) { # modification of image() from sp package # allows use of zlim if(missing(zlim)) zlim <- range(x@data[,1], na.rm=TRUE) plot(as(x, "Spatial"), xlim = NULL, ylim = NULL, axes = FALSE, asp = NA) image(as.image.SpatialGridDataFrame(x), add = TRUE, col = col, zlim = zlim) } # change prompt options(prompt="R> ") ################################################### ### chunk number 2: figlandsat ################################################### # display landsat data july1.ar <- radiocorr(july1, Grescale=0.77569, Brescale=-6.20, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1997, method="apparentreflectance") july61.thermal <- thermalband(july61, band=61) july.cloud <- clouds(july1.ar, july61.thermal) png("methodsPDF-plotlandsat.png", width=5, height=7.5, units="in", res=300) par(mfrow=c(3,2)) par.mai.orig <- par()["mai"] par.mai2 <- par.mai.orig par.mai2[3] <- .4 par.mai2[4] <- .1 par(mai=c(.1, .1, .4, .1)) imageSGDF(july4, zlim=c(0, 181), col=gray(0:50/50)) title("a. July") imageSGDF(nov4, zlim=c(0, 181), col=gray(0:50/50)) title("b. November") par(mai = par.mai2) hist(july4@data[is.na(july.cloud@data[,1]),1], breaks=seq(1, 181, by=2), xlab="DN", main="c. July data without clouds") hist(nov4@data[,1], breaks=seq(1, 181, by=2), xlab="DN", main="d. November data") par(mai=c(.1, .1, .4, .1)) imageSGDF(july.cloud) title("e. July cloud mask") # nov is cloudfree so not shown image(dem, col=topo.colors(48)) title("f. Elevation") par(mai = par.mai.orig) dev.off() ################################################### ### chunk number 3: georefshow ################################################### # using November as the reference image to calculate the needed shift july.shift <- georef(nov3, july3, maxdist=50) # adjust the band 1 image from July using the shift coefficients july1.corr <- geoshift(july1, padx=10, pady=10, july.shift$shiftx, july.shift$shifty) ################################################### ### chunk number 4: toarreflshow ################################################### # top of atmosphere reflectance example july4.ar <- radiocorr(july4, Grescale=0.63725, Brescale=-5.1, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1039, method="apparentreflectance") ################################################### ### chunk number 5: toarreflhide ################################################### # top of atmosphere reflectance - rest of images july1.ar <- radiocorr(july1, Grescale=0.77569, Brescale=-6.20, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1997, method="apparentreflectance") july2.ar <- radiocorr(july2, Grescale=0.79569, Brescale=-6.40, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1812, method="apparentreflectance") july3.ar <- radiocorr(july3, Grescale=0.61922, Brescale=-5.00, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1533, method="apparentreflectance") july5.ar <- radiocorr(july5, Grescale=0.12573, Brescale=-1.00, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=230.8, method="apparentreflectance") july7.ar <- radiocorr(july7, Grescale=0.04373, Brescale=-0.35, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=84.90, method="apparentreflectance") nov1.ar <- radiocorr(nov1, Grescale=0.77569, Brescale=-6.20, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=1997, method="apparentreflectance") nov2.ar <- radiocorr(nov2, Grescale=0.79569, Brescale=-6.40, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=1812, method="apparentreflectance") nov3.ar <- radiocorr(nov3, Grescale=0.61922, Brescale=-5.00, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=1533, method="apparentreflectance") nov4.ar <- radiocorr(nov4, Grescale=0.63725, Brescale=-5.1, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=1039, method="apparentreflectance") nov5.ar <- radiocorr(nov5, Grescale=0.12573, Brescale=-1.00, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=230.8, method="apparentreflectance") nov7.ar <- radiocorr(nov7, Grescale=0.04373, Brescale=-0.35, sunelev=26.2, edist=ESdist("2002-11-25"), Esun=84.90, method="apparentreflectance") ################################################### ### chunk number 6: thermalshow ################################################### # thermal band example july61.thermal <- thermalband(july61, band=61) ################################################### ### chunk number 7: thermalhide ################################################### # top of atmosphere reflectance example # uses default parameters july62.thermal <- thermalband(july62, band=62) nov61.thermal <- thermalband(nov61, band=61) nov62.thermal <- thermalband(nov62, band=62) ################################################### ### chunk number 8: cloudshow ################################################### july.cloud <- clouds(july1.ar, july61.thermal) nov.cloud <- clouds(nov1.ar, nov61.thermal) ################################################### ### chunk number 9: showrelnorm ################################################### # correct july to november, using band 4 as an example # using DN, removing the clouds in the July data july4.rncoef <- relnorm(nov4, july4, mask=july.cloud, nperm=0) july4.rn <- july4.rncoef[["newimage"]] ################################################### ### chunk number 10: showhistmatch ################################################### # correct july to november, using band 4 as an example # using DN, removing the clouds in the July data july4.hmcoef <- histmatch(nov4, july4, mask=july.cloud) july4.hm <- july4.hmcoef[["newimage"]] ################################################### ### chunk number 11: figstat ################################################### # transformed data png("methodsPDF-plotstat.png", width=5, height=5.33, units="in", res=300) par(mfrow=c(2,2)) par(cex = 0.66) par.mai.orig <- par()["mai"] par.mai2 <- par.mai.orig par.mai2[3] <- .4 par.mai2[4] <- .1 par(mai=c(.1, .1, .4, .1)) imageSGDF(july4.hm, zlim=c(10, 120), col=gray(0:50/50)) title("a. Histogram matching") imageSGDF(july4.rn, zlim=c(30, 71), col=gray(0:50/50)) title("b. Relative normalization") par(mai = par.mai2) hist(july4.hm@data[is.na(july.cloud@data[,1]),1], breaks=seq(10, 120, by=2), xlab="DN", main="c. Histogram matching") hist(july4.rn@data[is.na(july.cloud@data[,1]),1], breaks=seq(30, 71, by=1), xlab="DN", main="d. Relative normalization") dev.off() ################################################### ### chunk number 12: showpif1 ################################################### # identify PIFs for nov, the clearer image nov.PIF <- PIF(nov3, nov4, nov7, level=.90) # use major axis regression to correct band 4 for July july4.pifcorr <- lmodel2(nov4@data[nov.PIF@data[,1] == 1, 1] ~ july4@data[nov.PIF@data[,1] == 1, 1]) july4.pifcorr <- unlist(july4.pifcorr[["regression.results"]][2, 2:3]) ################################################### ### chunk number 13: showpif2 ################################################### # keep SpatialGridDataFrame format july4.pif <- july4 july4.pif@data[,1] <- july4@data[,1] * july4.pifcorr[2] + july4.pifcorr[1] # get rid of cloud areas july4.pif@data[!is.na(july.cloud@data[,1]), 1] <- NA ################################################### ### chunk number 14: showrcs ################################################### # identify RCSs for nov, the clearer image nov.tasscap <- tasscap("nov", 7) nov.RCS <- RCS(nov.tasscap) # use major axis regression to correct band 4 for July july4.rcscorr <- lmodel2(nov4@data[nov.RCS@data[,1] == 1, 1] ~ july4@data[nov.RCS@data[,1] == 1, 1]) july4.rcscorr <- unlist(july4.rcscorr[["regression.results"]][2, 2:3]) # keep SpatialGridDataFrame format july4.rcs <- july4 july4.rcs@data[,1] <- july4@data[,1] * july4.rcscorr[2] + july4.rcscorr[1] # get rid of cloud areas july4.rcs@data[!is.na(july.cloud@data[,1]), 1] <- NA ################################################### ### chunk number 15: figpif ################################################### # transformed data png("methodsPDF-plotpif.png", width=5, height=5.33, units="in", res=300) par(mfrow=c(2,2)) par(cex = 0.66) par.mai.orig <- par()["mai"] par.mai2 <- par.mai.orig par.mai2[3] <- .4 par.mai2[4] <- .1 par(mai=c(.1, .1, .4, .1)) imageSGDF(july4.pif, zlim=c(49, 52), col=gray(0:50/50)) title("a. Pseudo-invariant features") imageSGDF(july4.rcs, zlim=c(57, 74), col=gray(0:50/50)) title("b. Radiometric control set") par(mai = par.mai.orig) hist(july4.pif@data[,1], breaks=seq(49, 52, by=.1), xlab="DN", main="c. Pseudo-invariant features") hist(july4.rcs@data[,1], breaks=seq(57, 74, by=.5), xlab="DN", main="d. Radiometric control set") dev.off() ################################################### ### chunk number 16: showDOS1 ################################################### # identify the lowest DN with at least 1000 pixels # for use as Starting haze value SHV <- table(july1@data[,1]) SHV <- min(as.numeric(names(SHV)[SHV > 1000])) SHV ################################################### ### chunk number 17: showDOS2 ################################################### # use DOS() to calculate the SHV values for remaining bands # DNfinal.mean are the SHV values as calculated by Chavez 1989 july.DOS <- DOS(sat=7, SHV=SHV, SHV.band=1, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"))[["DNfinal.mean"]] july.DOS ################################################### ### chunk number 18: showDOScorr ################################################### # select the set of SHV appropriate for atmospheric conditions july.DOS <- july.DOS[,2] # use the appropriate SHV to correct bands 2 and 4 july2.DOSrefl <- radiocorr(july2, Grescale=0.79569, Brescale=-6.40, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1812, Lhaze=july.DOS[2], method="DOS") july4.DOSrefl <- radiocorr(july4, Grescale=0.63725, Brescale=-5.10, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1039, Lhaze=july.DOS[4], method="DOS") # get rid of cloud areas july4.DOSrefl@data[!is.na(july.cloud@data[,1]), 1] <- NA ################################################### ### chunk number 19: showCOSTZ ################################################### july4.COSTZrefl <- radiocorr(july4, Grescale=0.63725, Brescale=-5.10, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1039, Lhaze=july.DOS[4], method="COSTZ") # get rid of cloud areas july4.COSTZrefl@data[!is.na(july.cloud@data[,1]), 1] <- NA ################################################### ### chunk number 20: showDOS4 ################################################### july4.DOS4refl <- radiocorr(july4, Grescale=0.63725, Brescale=-5.10, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1039, Lhaze=july.DOS[4], method="DOS4") # get rid of cloud areas july4.DOS4refl@data[!is.na(july.cloud@data[,1]), 1] <- NA ################################################### ### chunk number 21: figabs ################################################### # transformed data png("methodsPDF-plotabs.png", width=5, height=5.33, units="in", res=300) par(mfrow=c(2,2)) par(cex = 0.66) par.mai.orig <- par()["mai"] par(mai=c(.1, .1, .4, .1)) july4.arC <- july4.ar july4.arC@data[!is.na(july.cloud@data[,1]), 1] <- NA imageSGDF(july4.arC, zlim=c(0,1), col=gray(0:50/50)) title("a. At-sensor reflectance") imageSGDF(july4.DOSrefl, zlim=c(0,1), col=gray(0:50/50)) title("b. DOS") imageSGDF(july4.COSTZrefl, zlim=c(0,1), col=gray(0:50/50)) title("c. COSTZ") imageSGDF(july4.DOS4refl, zlim=c(0,1), col=gray(0:50/50)) title("d. Modified DOS") par(mai = par.mai.orig) dev.off() ################################################### ### chunk number 22: showtopo1 ################################################### dem.slopeasp <- slopeasp(dem) nov4.cosine <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="cosine") ################################################### ### chunk number 23: showtopo2 ################################################### nov4.cosimp <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="improvedcosine") ################################################### ### chunk number 24: showtopo8 ################################################### nov4.gamma <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="gamma") ################################################### ### chunk number 25: showtopo7 ################################################### nov4.SCS <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="SCS") ################################################### ### chunk number 26: showtopo3 ################################################### nov4.minnaert <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="minnaert") ################################################### ### chunk number 27: showtopo4 ################################################### nov4.min2 <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="minslope") ################################################### ### chunk number 28: figtopo1 ################################################### png("methodsPDF-plottopo1.png", width=5, height=7.5, units="in", res=300) par(mfrow=c(3,2)) par.mai.orig <- par()["mai"] par(mai=c(.1, .1, .4, .1)) imageSGDF(nov4.cosine, col=gray(0:50/50), zlim=c(0, .6)) title("a. Cosine correction") imageSGDF(nov4.cosimp, col=gray(0:50/50), zlim=c(0, .6)) title("b. Improved cosine correction") imageSGDF(nov4.gamma, col=gray(0:50/50), zlim=c(0, .6)) title("c. Gamma correction") imageSGDF(nov4.SCS, col=gray(0:50/50), zlim=c(0, .6)) title("d. SCS correction") imageSGDF(nov4.minnaert, col=gray(0:50/50), zlim=c(0, .6)) title("e. Minnaert correction") imageSGDF(nov4.min2, col=gray(0:50/50), zlim=c(0, .6)) title("f. Minnaert with slope correction") par(mai = par.mai.orig) dev.off() ################################################### ### chunk number 29: showtopo5 ################################################### nov4.ccor <- topocorr(nov4.ar, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=26.2, sunazimuth=159.5, method="ccorrection") ################################################### ### chunk number 30: showtopo6 ################################################### dem.smooth5 <- slopeasp(dem, smoothing=5) nov4.ccorsmooth <- topocorr(nov4.ar, dem.smooth5$slope, dem.smooth5$aspect, sunelev=26.2, sunazimuth=159.5, method="ccorrection") ################################################### ### chunk number 31: figtopo2 ################################################### png("methodsPDF-plottopo2.png", width=5, height=5.33, units="in", res=300) par(mfrow=c(2,2)) par(cex = 0.66) par(mai=c(.1, .1, .4, .1)) par.mai.orig <- par()["mai"] imageSGDF(nov4.ccor, col=gray(0:50/50), zlim=c(0, .6)) title("a. C-correction") imageSGDF(nov4.ccorsmooth, col=gray(0:50/50), zlim=c(0, .6)) title("b. C-correction with smoothing") imageSGDF(dem.slopeasp$slope, col=terrain.colors(24), zlim=c(0, 32)) title("c. Slope") imageSGDF(dem.smooth5$slope, col=terrain.colors(24), zlim=c(0, 32)) title("d. Slope with a smoothing factor of 5") par(mai = par.mai.orig) dev.off()