R

ランドサット8号のデータ(反射率のプロット)

rLandsat8 , RStoolbox , raster パッケージ

場所は鳥取市周辺

  • 画像上の4点の反射率をバンド別にプロット

地点の座標を見つけるのが面倒。こういうことするのはGUIのソフトがいい。

code

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
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)
#
#rLandsat8 packageの関数の一部訂正と付け加え
#反射率（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")
#鳥取市 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)
#リスト l を作成
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)
#
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"))
#反射率をプロットする地点
pos1<-c(423000,3930000)
pos2<-c(423000,3927400)
pos3<-c(429650,3933400)
pos4<-c(424240,3932130)
#画像と反射率をプロットする地点
#png("ref_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()
#
#仮に複数のデータがあっても最小値を取り出す。p1@data@min
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))
#png("ref_tottori02.png",width=800,height=800)
plot(1:length(p1),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)
#dev.off()