R

ランドサット8号のデータ(大気補正)

landsat , RStoolbox , raster パッケージ

場所は鳥取市周辺

  • 大気補正した後画像上の4点の反射率をバンド別にプロット
  1. landsat パッケージのradiocorrを元にした
  2. esunはここの数値を使った。
  3. radiance_darkは画像の最も暗いピクセルの数値

（参考）
Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects?

今回作成した関数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#method:DOS1 , DOS2 , DOS4
#ランドサット8号の場合
#DOS2:TAUz=suntheta "aerosol", "blue", "green", "red", "nir"
#     TAUz=1        "swir1", "swir2"
radiocorr<-function(landsat8, band,method = "DOS1") {
    bandnames <- c("aerosol", "blue", "green", "red", "nir", 
        "swir1", "swir2")
  idx <- seq_along(bandnames)[sapply(bandnames, function(x) band %in% x)]
  ml <- as.numeric(landsat8$metadata[[paste0("radiance_mult_band_",idx)]])
  al <- as.numeric(landsat8$metadata[[paste0("radiance_add_band_",idx)]])
  TOArad <- landsat8$band[[band]] * ml + al
sunElev<-as.numeric(landsat8$metadata$sun_elevation)
suntheta <- cos((90 - sunElev) * pi / 180)
satzenith=as.numeric(landsat8$metadata$roll_angle)
    satphi <- cos(satzenith * pi/180)
d<-as.numeric(landsat8$metadata$earth_sun_distance)
#esun
#http://sourceforge.net/p/saga-gis/code-0/2166/tree//trunk/saga-gis/src/modules/imagery/imagery_tools/landsat_toar_core.cpp#l1016
#{ 2026.8, 2066.8, 1892.5, 1602.8, 972.6, 245.0, 79.7, 1805.5, 399.7, 0., 0. }
#esun<-NULL
# 大気補正 haze:1%
  if (band=="aerosol") {esun=2026.8}
  if (band=="blue") {esun=2066.8}
  if (band=="green") {esun=1892.5} 
  if (band=="red") {esun=1602.8}
  if (band=="nir") {esun=972.6}    
  if (band=="swir1") {esun=245.0} 
  if (band=="swir2") {esun=79.7}
    if (method =="DOS1") {
        TAUz <- 1
        TAUv <- 1
        Edown <- 0
        Lp <- TOArad@data@min-0.01*(esun*suntheta)/(d^2*pi)
    }
    if (method =="DOS2") {
        TAUz <- suntheta
        TAUv <- 1
        Edown <- 0
        Lp <- TOArad@data@min-0.01*(esun*suntheta)/(d^2*pi)
    }
    else if (method =="DOS4") {
        TAUv <- 1 ; TAUz <- 1
        taudiff <- 1
        tau <- 9999
        Edown <- 0
        while (abs(taudiff) > 1e-07) {
            taudiff <- tau
#            Eo <- esun/d^2
            Lp <- TOArad@data@min-0.01 *(esun * suntheta *TAUz+Edown)*TAUv/(d^2*pi)   
            taustep <- 1 - (4 * pi * Lp)/(esun * suntheta/d^2)
            while (taustep < 0) {
                Lp <-TOArad@data@min-0.01*(esun*suntheta*TAUz+Edown)*TAUv/(d^2*pi)
                taustep <- 1 - (4 * pi * Lp)/(esun * suntheta/d^2)
            }
            tau <- -1 * suntheta * log(1 - (4 * pi * Lp)/(esun *suntheta/d^2))
            TAUv <- exp(-1 * tau/satphi)
            TAUz <- exp(-1 * tau/suntheta)
            Edown <- pi * Lp
            taudiff <- taudiff - tau
        }
}
    results <- (pi * d^2 *(TOArad-Lp))/(TAUv*(esun*suntheta*TAUz+Edown))
    return(results)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
library(raster)
library(colorRamps)
library(RColorBrewer)
#library(devtools)
#install_url("https://github.com/Terradue/rLandsat8/releases/download/v0.1-SNAPSHOT/rLandsat8_0.1.0.tar.gz")
library(rLandsat8)
#
#反射率（Reflectance）
#太陽天頂角を考慮する（変更後）
ToTOAReflectance <- function(landsat8, band) {
    bandnames <- c("aerosol", "blue", "green", "red", "nir", 
        "swir1", "swir2", "panchromatic", "cirrus", "tirs1", 
        "tirs2")
    idx <- seq_along(bandnames)[sapply(bandnames, function(x) band %in% 
        x)]
    ml <- as.numeric(landsat8$metadata[[paste0("reflectance_mult_band_", 
        idx)]])
    al <- as.numeric(landsat8$metadata[[paste0("reflectance_add_band_", 
        idx)]])
 #   TOAref <- landsat8$band[[band]] * ml + al
    TOAref <- (landsat8$band[[band]] * ml + al)/sin(as.numeric(landsat8$metadata$sun_elevation) * pi /180)
    return(TOAref)
}
#
#改良型正規化水指数（MNDWI)
#greenがgreebになっていたのを直した。
#
ToMNDWI <- function(landsat8) {
  # MNDWI=(Green - SWIR1)/(Green + SWIR1) 
  green <- ToTOAReflectance(landsat8, "green")
  swir1 <- ToTOAReflectance(landsat8, "swir1")
# mndwi <- (greeb - swir1)/(green + swir1)
  mndwi <- (green - swir1) / (green + swir1)
  return(mndwi)
}
#
#付け加えた関数（landsat8のメタファイル(ファイル名の後ろが"_MTL.txt")を読み込む）
#
read.meta<- function(product){
meta.file <- paste0(product, "/", product, "_MTL.txt")
    textLines <- readLines(meta.file)
    met <- read.table(text=textLines[grep("=", textLines)],header = FALSE, sep = "=", strip.white = TRUE, 
    stringsAsFactors = FALSE,row.names = NULL, col.names = c("name", "value"))
    met <- met[!met$name == "GROUP", ]
    met <- met[!met$name == "END_GROUP", ]
    rownames(met) <- tolower(met[, "name"])
    met[, "name"] <- NULL
    met <- as.list(as.data.frame(t(met), stringsAsFactors = FALSE))
return(met)       
}
#
##環境によってはRStoolboxパッケージがインストールできないのでpanSharpenのみ取り出し。
#使う方法は　method = "brovey"だけにした。
panSharpen <- function(img, pan, r, g, b) {
	stopifnot(inherits(img, "Raster") & inherits(pan, "Raster"))
        layernames <- names(img)[c(r,g,b)] 
        ordr <- match(layernames, names(img))
        layernames <- layernames[order(ordr)]
        msi    <- resample(img[[layernames]], pan, method = "ngb")
        mult   <- pan / sum(msi)
        panimg <- msi * mult
    names(panimg) <- paste0(layernames, "_pan")
    panimg
}
#
#
#setwd("/home/user/landsat8")
setwd("/media/user/MUSIC")
#setwd("e:/")
#鳥取市 path:111 row:035
product  <- "LC81110352015214LGN00"
Crop <- c(420000,430000,3926000,3934600)
#fname0<-"tottori"
#メタデータの読み込み
meta<-read.meta(product)
date<-meta$date_acquired
#保存するファイル名,作成する画像のファイル名の基礎
fname<-paste0(fname0,substr(date,1,4),substr(date,6,7),substr(date,9,10))
#解像度=30m
OLI<-stack(c(aerosol=paste0(product,"/",product, "_B1.TIF"),blue=paste0(product,"/",product, "_B2.TIF"),
green=paste0(product,"/",product, "_B3.TIF"),red=paste0(product,"/",product, "_B4.TIF"),
nir=paste0(product,"/",product, "_B5.TIF"),swir1=paste0(product,"/",product, "_B6.TIF"),
swir2=paste0(product,"/",product, "_B7.TIF"),cirrus=paste0(product,"/",product, "_B9.TIF")))
#解像度=15m
PAN <-stack(paste0(product,"/",product, "_B8.TIF"))
names(PAN)<-"panchromatic"
#解像度=100m
TIRS<-stack(c(tirs1=paste0(product,"/",product, "_B10.TIF"),tirs2=paste0(product,"/",product, "_B11.TIF")))
#必要な部分の切り出し
subOLI <- crop(OLI,Crop)
subPAN <- crop(PAN,Crop)
subTIRS <- crop(TIRS,Crop)
#
bands=c(unstack(subOLI)[1:7],unstack(subPAN),unstack(subOLI)[8],unstack(subTIRS))
l<-list(metadata = meta, band = bands)
names(l$band)<-c("aerosol","blue", "green", "red","nir","swir1","swir2","panchromatic","cirrus","tirs1","tirs2")
#
RGBpan <- panSharpen(stack(l$band[2:4]),l$band$panchromatic, r = 3, g = 2, b = 1)
#
pos1<-c(423000,3930000)
pos2<-c(423000,3927400)
pos3<-c(429650,3933400)
pos4<-c(424240,3932130)
#
#png("refcor_tottori01.png",width=800,height=800)
plotRGB(RGBpan,r = 3, g = 2, b = 1,stretch="lin")
points(pos1[1],pos1[2],col="red",pch=16,cex=1.5)
points(pos2[1],pos2[2],col="blue",pch=16,cex=1.5)
points(pos3[1],pos3[2],col="green",pch=16,cex=1.5)
points(pos4[1],pos4[2],col="orange",pch=16,cex=1.5)
#dev.off()
#
#大気補正なし
l_ref<-stack(ToTOAReflectance(landsat8=l, band="aerosol"),
	ToTOAReflectance(landsat8=l, band="blue"),
	ToTOAReflectance(landsat8=l, band="green"),
	ToTOAReflectance(landsat8=l, band="red"),
	ToTOAReflectance(landsat8=l, band="nir"),
	ToTOAReflectance(landsat8=l, band="swir1"),
	ToTOAReflectance(landsat8=l, band="swir2"))
#大気補正あり
#DOS1
l_ref_cor<-stack(radiocorr(l,band="aerosol",method = "DOS1"),
        radiocorr(l,band="blue",method = "DOS1"),
        radiocorr(l,band="green",method = "DOS1"),
        radiocorr(l,band="red",method = "DOS1"),
        radiocorr(l,band="nir",method = "DOS1"),
        radiocorr(l,band="swir1",method = "DOS1"),
        radiocorr(l,band="swir2",method = "DOS1"))
#DOS2
l_ref_dos2<-stack(radiocorr(l,band="aerosol",method = "DOS2"),
        radiocorr(l,band="blue",method = "DOS2"),
        radiocorr(l,band="green",method = "DOS2"),
        radiocorr(l,band="red",method = "DOS2"),
        radiocorr(l,band="nir",method = "DOS2"),
        radiocorr(l,band="swir1",method = "DOS1"),
        radiocorr(l,band="swir2",method = "DOS1"))
#DOS4
l_ref_dos4<-stack(radiocorr(l,band="aerosol",method = "DOS4"),
        radiocorr(l,band="blue",method = "DOS4"),
        radiocorr(l,band="green",method = "DOS4"),
        radiocorr(l,band="red",method = "DOS4"),
        radiocorr(l,band="nir",method = "DOS4"),
        radiocorr(l,band="swir1",method = "DOS4"),
        radiocorr(l,band="swir2",method = "DOS4"))
#
p1 <- crop(l_ref,c(pos1[1],pos1[1]+1,pos1[2],pos1[2]+1))
p2 <- crop(l_ref,c(pos2[1],pos2[1]+1,pos2[2],pos2[2]+1))
p3 <- crop(l_ref,c(pos3[1],pos3[1]+1,pos3[2],pos3[2]+1))
p4 <- crop(l_ref,c(pos4[1],pos4[1]+1,pos4[2],pos4[2]+1))
#大気補正あり(DOS1)
p5 <- crop(l_ref_cor,c(pos1[1],pos1[1]+1,pos1[2],pos1[2]+1))
p6 <- crop(l_ref_cor,c(pos2[1],pos2[1]+1,pos2[2],pos2[2]+1))
p7 <- crop(l_ref_cor,c(pos3[1],pos3[1]+1,pos3[2],pos3[2]+1))
p8 <- crop(l_ref_cor,c(pos4[1],pos4[1]+1,pos4[2],pos4[2]+1))
#大気補正あり(DOS4)
p9 <- crop(l_ref_dos4,c(pos1[1],pos1[1]+1,pos1[2],pos1[2]+1))
p10 <- crop(l_ref_dos4,c(pos2[1],pos2[1]+1,pos2[2],pos2[2]+1))
p11 <- crop(l_ref_dos4,c(pos3[1],pos3[1]+1,pos3[2],pos3[2]+1))
p12 <- crop(l_ref_dos4,c(pos4[1],pos4[1]+1,pos4[2],pos4[2]+1))
#大気補正あり(DOS2)
p13 <- crop(l_ref_dos2,c(pos1[1],pos1[1]+1,pos1[2],pos1[2]+1))
p14 <- crop(l_ref_dos2,c(pos2[1],pos2[1]+1,pos2[2],pos2[2]+1))
p15 <- crop(l_ref_dos2,c(pos3[1],pos3[1]+1,pos3[2],pos3[2]+1))
p16 <- crop(l_ref_dos2,c(pos4[1],pos4[1]+1,pos4[2],pos4[2]+1))
#png("refcor_tottori02.png",width=1000,height=1000)
par(mfrow=c(2,2))
plot(p1@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
title("場所によるバンド別反射率の違い(大気補正なし)")
axis(1, at=1:length(p1), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
lines(p2@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue")
lines(p3@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green")
lines(p4@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange")
legend("topleft",legend=c("池","森林","砂浜","滑走路"),col=c("red","blue","green","orange"),pch=16,cex=1,lty=1)
#
plot(p5@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
title("場所によるバンド別反射率の違い(大気補正あり:DOS1)")
axis(1, at=1:length(p5), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
lines(p6@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue")
lines(p7@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green")
lines(p8@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange")
legend("topleft",legend=c("池","森林","砂浜","滑走路"),col=c("red","blue","green","orange"),pch=16,cex=1,lty=1)
#
plot(p13@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
title("場所によるバンド別反射率の違い(大気補正あり:DOS2)")
axis(1, at=1:length(p1), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
lines(p14@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue")
lines(p15@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green")
lines(p16@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange")
legend("topleft",legend=c("池","森林","砂浜","滑走路"),col=c("red","blue","green","orange"),pch=16,cex=1,lty=1)
#
plot(p9@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
title("場所によるバンド別反射率の違い(大気補正あり:DOS4)")
axis(1, at=1:length(p9), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
lines(p10@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue")
lines(p11@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green")
lines(p12@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange")
legend("topleft",legend=c("池","森林","砂浜","滑走路"),col=c("red","blue","green","orange"),pch=16,cex=1,lty=1)
par(mfrow=c(1,1))
#dev.off()
#png("refcor_tottori03.png",width=1000,height=1000)
par(mfrow=c(2,2))
plot(p1@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
axis(1, at=1:length(p1), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
title("大気補正によるバンド別反射率の違い(池)")
lines(p5@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue",lty=2)
lines(p9@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green",lty=3)
lines(p13@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange",lty=4)
legend("topleft",legend=c("大気補正なし","ＤＯＳ１","ＤＯＳ２","ＤＯＳ４"),col=c("red","blue","orange","green"),pch=16,cex=1,lty=c(1,2,4,3))
#
plot(p2@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
axis(1, at=1:length(p2), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
title("大気補正によるバンド別反射率の違い(森林)")
lines(p6@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue",lty=2)
lines(p10@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green",lty=3)
lines(p14@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange",lty=4)
legend("topleft",legend=c("大気補正なし","ＤＯＳ１","ＤＯＳ２","ＤＯＳ４"),col=c("red","blue","orange","green"),pch=16,cex=1,lty=c(1,2,4,3))
#
plot(p3@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
axis(1, at=1:length(p3), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
title("大気補正によるバンド別反射率の違い(砂浜)")
lines(p7@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue",lty=2)
lines(p11@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green",lty=3)
lines(p15@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange",lty=4)
legend("topleft",legend=c("大気補正なし","ＤＯＳ１","ＤＯＳ２","ＤＯＳ４"),col=c("red","blue","orange","green"),pch=16,cex=1,lty=c(1,2,4,3))
#
plot(p4@data@min,type="o",pch=16,cex=1.5,lwd=2,col="red",xaxt = "n",
xlab="Band",ylab="Reflectance",las=1,ylim=c(0,0.3),bty = "l")
axis(1, at=1:length(p4), labels=c("aerosol\n0.43-0.45μm","blue\n0.45-0.51μm", "green\n0.53-0.59μm",
 "red\n0.64-0.67μm","nir\n0.85-0.88μm","swir1\n1.57-1.65μm","swir2\n2.11-2.29μm"),tck=0.01)
title("大気補正によるバンド別反射率の違い(滑走路)")
lines(p8@data@min,type="o",pch=16,cex=1.5,lwd=2,col="blue",lty=2)
lines(p12@data@min,type="o",pch=16,cex=1.5,lwd=2,col="green",lty=3)
lines(p16@data@min,type="o",pch=16,cex=1.5,lwd=2,col="orange",lty=4)
legend("topleft",legend=c("大気補正なし","ＤＯＳ１","ＤＯＳ２","ＤＯＳ４"),col=c("red","blue","orange","green"),pch=16,cex=1,lty=c(1,2,4,3))
par(mfrow=c(1,1))
#dev.off()