R

ランドサット8号から見た西之島(code)

raster , rasterVis パッケージ

データはLandBrowserより

True画像と温度を推定する画像(band10とband7)

パッケージと関数を読み込む
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
library(raster)
library(rasterVis)
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)
#source("rLandsat8fix.R") #変更した関数付け加えた関数をrLandsat8fix.Rというファイル名で作業フォルダに保存した場合。
#付け加えた関数（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
}
#反射率（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)
}
#反射率（Reflectance）大気補正：method DOS1,DOS2,DOS4
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)
}
15個のデータを繰り返し処理
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
setwd("/media/user/MUSIC/nishinoshima")
dir<-c("LC81050412013310LGN00","LC81050412013358LGN00","LC81050412014089LGN00","LC81050412014185LGN00",
	"LC81050412014281LGN00","LC81050412014297LGN00","LC81050412014329LGN00","LC81050412015044LGN00",
	"LC81050412015060LGN00","LC81050412015092LGN00","LC81050412015156LGN00","LC81050412015172LGN00",
	"LC81050412015188LGN00","LC81050412015204LGN00","LC81050412015236LGN00")
fname0<-"nishinoshima"
#15個のデータを繰り返し処理
for (i in 1:15){
#i=5
#西之島 path:105 row:041
product  <- dir[i]
Crop<-c(485700,491000,3011000,3016000)
#メタデータの読み込み
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")
#save(l,file=paste0(fname,"_1.RData"))
## Brovey pan sharpening
RGBpan <- panSharpen(stack(l$band[2:4]),l$band$panchromatic,r=3,g=2,b=1)
## Plot 
png(paste0(fname,"_1.png"),width=1200,height=1000)
par(mfrow=c(1,2))
plotRGB(RGBpan,r = 3, g = 2, b = 1,stretch = "hist") 
text((Crop[1]+Crop[2])/2,Crop[4],"Pan-sharpened(method：brovey stretch:hist)",pos=3, cex=1.5)
plotRGB(RGBpan,r = 3, g = 2, b = 1,stretch = "lin") 
text((Crop[1]+Crop[2])/2,Crop[4],"Pan-sharpened(method：brovey stretch:lin)",pos=3, cex=1.5)
par(mfrow=c(1,1))
dev.off()
#setwd("/home/user/landsat8")
#load(paste0(fname,"_1.RData")) #保存したデータ名
#スペクトル放射輝度から温度(btemp: 絶対温度(K))
btemp <- ToAtSatelliteBrightnessTemperature(l, band="tirs1")
# 絶対温度を摂氏に変換
SST <- btemp - 273.15  
temperature<-radiocorr(landsat8=l, band="swir2",method="DOS4")
#輝度温度の色分けは10～80で固定
png(paste0(fname,"_2.png"),width=1200,height=1000)
par(mfrow=c(1,2))
plot(l$band$panchromatic, col = gray.colors(30),legend=FALSE, axes=FALSE,box=FALSE)
#plot(SST,col=blue2green2red(30),legend=T,axes=FALSE,box=FALSE,alpha=0.5,add=T)
breakpoints <- seq(10,80,2)
plot(SST,col=blue2green2red(length(breakpoints)),legend=F,axes=FALSE,box=FALSE,alpha=0.5,add=T,breaks=breakpoints)
plot(SST, legend.only=TRUE,col=blue2green2red(length(breakpoints)),
     legend.width=1, legend.shrink=0.5,,breaks=breakpoints,
     axis.args=list(at=seq(10,80,5),
                    labels=seq(10,80,5), 
                    cex.axis=1))
text((Crop[1]+Crop[2])/2,Crop[4],"輝度温度（℃）",pos=3, cex=1.5)
plot(l$band$panchromatic, col = gray.colors(30),legend=FALSE, axes=FALSE,box=FALSE)
plot(temperature,col=blue2green2red(30),legend=F,axes=FALSE,box=FALSE,alpha=0.5,add=T)
r.range <- c(minValue(temperature), maxValue(temperature))
plot(temperature, legend.only=TRUE, col=blue2green2red(30),
     legend.width=1, legend.shrink=0.5,
     axis.args=list(at=c(r.range[1], r.range[2]),
                    labels=c("L","H"), 
                    cex.axis=1),
     legend.args=list(text='temperature', side=4, font=1, line=0.75, cex=1))
text((Crop[1]+Crop[2])/2,Crop[4],"BAND7(swir2)",pos=3, cex=1.5)
par(mfrow=c(1,1))
dev.off()
}
輝度温度が８０℃を越えるデータをrasterVis::levelplotで視覚化
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
#i=2
#i=9
#i=13
i=14
#西之島 path:105 row:041
product  <- dir[i]
Crop<-c(485700,491000,3011000,3016000)
#メタデータの読み込み
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")
#スペクトル放射輝度から温度(btemp: 絶対温度(K))
btemp <- ToAtSatelliteBrightnessTemperature(l, band="tirs1")
SST <- btemp - 273.15  # 絶対温度を摂氏に変換
png(paste0(fname,"_3.png"),width=800,height=800)
levelplot(SST,contour=F,margin=T,at=seq(10,100,1),par.settings=BuRdTheme)
dev.off()