R

ランドサット8号のデータ3(クラスタリング)

rLandsat8 , raster , sp , colorRamps パッケージ

「ランドサット8号のデータ」のまとめ

ランドサット８号のデータを入手
LandBrowser
データ検索＆ダウンロード

Path/RowとLatitude/Longitudeの変換
WRS-2 Path/Row to Latitude/Longitude Converter

緯度経度座標をＵＴＭ座標に変換する

Rを使って緯度経度座標をＵＴＭ座標に変換する

1
2
3
4
5
6
library(rgdal)
cord.dec = SpatialPoints(cbind(130.59,31.545), proj4string=CRS("+proj=longlat"))
spTransform(cord.dec, CRS("+proj=utm +zone=52"))
##
cord.dec = SpatialPoints(cbind(130.73,31.63), proj4string=CRS("+proj=longlat"))
spTransform(cord.dec, CRS("+proj=utm +zone=52"))

ダウンロードしたデータを解凍して作業フォルダに置く（ここでは /home/user/landsat8 とする）

1
2
3
4
5
6
7
8
setwd("/home/user/landsat8")
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)
#source("rLandsat8fix.R") #変更した関数付け加えた関数をrLandsat8fix.Rというファイル名で作業フォルダに保存。

rLandsat8パッケージの関数を少し変更

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
#反射率（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)       
}

最初にデータの一部を切り出してからデータ処理をする。
記事「ランドサット8号のデータ1」のやり方より速い。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
product  <- "LC81110352015214LGN00"
#メタデータを読み込む
meta<-read.meta(product)
#例えば、utm_zoneは as.numeric(meta$utm_zone)
#鳥取市の切り出し場所
#「緯度経度座標をＵＴＭ座標に変換」を参考にする
Crop <- c(420000,430000,3926000,3934600)
#解像度=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)

可視光の地図

1
2
3
4
5
6
#png("tottori20150802_1.png",width=1200,height=1000)
par(mfrow=c(1,2))
plotRGB(subOLI,4,3,2)
plotRGB(subOLI,4,3,2,stretch="lin")
par(mfrow=c(1,1))
#dev.off()

rLandsat8パッケージの関数の関数が使えるようにする

1
2
3
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")
放射輝度（Radiance）
1
2
3
#radiance.blue <- ToTOARadiance(landsat8=l, band="blue")
#radiance.green <- ToTOARadiance(landsat8=l, band="green")
#radiance.red <- ToTOARadiance(landsat8=l, band="red")
反射率（Reflectance）

太陽天頂角を考慮しない(変更前)

class : RasterLayer
dimensions : 287, 333, 95571 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 420015, 430005, 3925995, 3934605 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : blue
values : 0.09032, 0.3377 (min, max)

太陽天頂角を考慮する（変更後）

1
2
3
#sin(as.numeric(l$metadata$sun_elevation) * pi /180) 
#0.8912265  <-この数で割る。
ToTOAReflectance(landsat8=l, band="blue")

class : RasterLayer
dimensions : 287, 333, 95571 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 420015, 430005, 3925995, 3934605 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : blue
values : 0.1013435, 0.378916 (min, max)

スペクトル放射輝度から温度(btemp: 絶対温度(K))
1
2
3
4
5
6
7
btemp <- ToAtSatelliteBrightnessTemperature(l, band="tirs1")
SST <- btemp - 273.15  # 絶対温度を摂氏に変換
#google earthで扱うためにkmz(kml)形式で保存
require(rgdal)
newproj <- "+proj=longlat +ellps=WGS84"
SST_spdf <- projectRaster(SST,crs=newproj) 
KML(SST_spdf, file='tottori_temp',col=blue2green2red(30))

google earthに重ねると

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#正規化植生指数(NDVI)
ndvi<-ToNDVI(l)
#地表水分指標(LSWI)
lswi <- ToLSWI(l)
#改良型正規化水指数（MNDWI)
mndwi<-ToMNDWI(l)
#png("tottori20150802_2.png",width=1200,height=1000)
par(mfrow=c(2,2))
plot(btemp-273.15,col=blue2green2red(30)) 　#絶対温度を摂氏にするために273.15を引く
plot(ndvi)
plot(lswi)
plot(mndwi)
par(mfrow=c(1,1))
#dev.off()

  • 温度（左上）　正規化植生指数（右上）
  • 地表水分指標（左下）　改良型正規化水指数（右下）
クラスタリング

aerosol 、 blue、 green、 red、 nir、 swir1、swir2 を使う

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
aerosol<-ToTOAReflectance(landsat8=l,band="aerosol")
blue<-ToTOAReflectance(landsat8=l,band="blue")
green<-ToTOAReflectance(landsat8=l,band="green")
red<-ToTOAReflectance(landsat8=l,band="red")
nir<-ToTOAReflectance(landsat8=l,band="nir")
swir1<-ToTOAReflectance(landsat8=l,band="swir1")
swir2<-ToTOAReflectance(landsat8=l,band="swir2")
cl<- stack(aerosol,blue,green,red,nir,swir1,swir2) 
DF <- as.data.frame(cl) ;head(cl);tail(cl)
#クラスタリング cluster::clara関数
library(cluster)
E=clara(DF,8)
clusters <- raster(cl)   ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, E$clustering) 
clusters
#自己組織化マップによるクラスタリング som::som
library(som)
E.som<-som(DF, xdim=8, ydim=1)
#E.som$visual[,1]
clusters1 <- raster(cl)   ## create an empty raster with same extent than ICE
clusters1 <- setValues(clusters1, E.som$visual[,1]+1)  #番号を1から8にするために1加える
clusters1
#png("tottori20150802_3.png",width=1200,height=1000)
par(mfrow=c(1,2))
plot(clusters,col=brewer.pal(8, "Paired"))
plot(clusters1,col=brewer.pal(8, "Paired"))
par(mfrow=c(1,1))
#dev.off()

配色を変更

1
2
3
4
5
6
#png("tottori20150802_4.png",width=1200,height=1000)
par(mfrow=c(1,2))
plot(clusters,col=c("#A6CEE3","#E31A1C","#FB9A99","#FDBF6F" ,"#FF7F00","#B2DF8A","#33A02C","#1F78B4" ))
plot(clusters1,col=c("#A6CEE3", "#1F78B4","#E31A1C","#FDBF6F","#FB9A99","#FF7F00", "#B2DF8A","#33A02C"))
par(mfrow=c(1,1))
#dev.off()

自己組織化マップによるクラスタリング の結果をgoogle earthで扱うためにkmz(kml)形式で保存 

1
2
3
newproj <- "+proj=longlat +ellps=WGS84"
clusters1_spdf <- projectRaster(clusters1,crs=newproj) 
KML(clusters1_spdf, file='tottori_cls',col=c("#A6CEE3", "#1F78B4","#E31A1C","#FDBF6F","#FB9A99","#FF7F00", "#B2DF8A","#33A02C"))

google earthに重ねると