R

ランドサット8号のデータ

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

rLandsat8 パッケージを使ってみた（バージョンは0.1.0）

(参考)
ランドサット8号の日本上空からのデータを即時公開

LandBrowser
今回は、LC81110352015214LGN00.tar.bz2(2015-08-02のデータ)をダウンロードした。

(注意：ファイル名が～KUJ00.tar.bz2となっているファイル)
データ検索＆ダウンロード　AIST CS-W Client (デスクトップ用）
path:111-111,row:35-35で検索したデータ例えば
「LC81110352015214KUJ00.tar.bz2」 LGNではなくKUJとなっているデータは B10.TIF B11.TIF がない。
読み込めなかった。
Acquisition StationをLGNにしても No data to display !!

(rLandsat8 パッケージの使い方)
rLandsat8 - an R interface to Landsat8

インストール方法

1
2
#library(devtools)
#install_url("https://github.com/Terradue/rLandsat8/releases/download/v0.1-SNAPSHOT/rLandsat8_0.1.0.tar.gz")

(準備)
適当なフォルダ　ここでは　/home/user/landsat8/LC81110352015214LGN00 の中に解凍したデータをすべて入れる 

1
2
3
4
5
6
7
8
9
10
11
12
13
library(rLandsat8)
setwd("/home/user/landsat8")
# LC81110352015214LGN00 はフォルダ名。
product  <- "LC81110352015214LGN00"
# LC81110352015214LGN00 内にある メタデータ、tif データを読み込む。
l <- ReadLandsat8(product)
#Access the product metadata
l$metadata #出力は省略
#path row
l$metadata$wrs_path #出力は省略
l$metadata$wrs_row #出力は省略
#バンド名
names(l$band)

“aerosol” “blue” “green” “red” “nir”
“swir1” “swir2” “panchromatic” “cirrus” “tirs1” “tirs2” 

1
2
3
4
5
6
7
8
9
10
#To access a band:
l$band$blue
l[[2]]$blue
#plot
#png("landsat8_001.png",width=1000,height=1000)
plot(l$band$tirs1)
#dev.off()
#plot(l[[2]]$blue)
#plot(l[[2]]$green)
#plot(l[[2]]$red)

メタデータから必要な数値を読み込むので計算式を作る必要がない

放射輝度
1
2
3
radiance.blue <- ToTOARadiance(landsat8=l, band="blue")
radiance.green <- ToTOARadiance(landsat8=l, band="green")
radiance.red <- ToTOARadiance(landsat8=l, band="red")
反射率
1
2
3
#no sun angle correction
reflectance.blue <- ToTOAReflectance(landsat8=l, band="blue")
reflectance.blue

class : RasterLayer
dimensions : 7901, 7761, 61319661 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 290385, 523215, 3870285, 4107315 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /tmp/R_raster_user/r_tmp_2015-08-16_105516_1892_04310.grd
names : LC81110352015214LGN00_B2
values : 0.06486, 0.925 (min, max) 

1
2
#using sun angle correction
reflectance.green <- ToTOAReflectance(landsat8=l, band="green", is.suncorrected = TRUE)
スペクトル放射輝度から輝度温度へ
1
2
3
4
5
btemp <- ToAtSatelliteBrightnessTemperature(l, band="tirs1")
#png("landsat8_002.png",width=1000,height=1000)
SST <- btemp - 273.15  # 絶対温度を摂氏に変換
plot(SST)
#dev.off()

NDVI (正規化植生指数)

1
2
3
4
5
#Normalized Difference Vegetation Index (NDVI)
ndvi <- ToNDVI(l)
#png("landsat8_003.png",width=1000,height=1000)
plot(ndvi)
#dev.off()

1
2
#Land Surface Water Index (LSWI)
LSWI <- ToLSWI(l)

動かなかったコード

1
2
3
4
5
6
7
8
9
#using sun angle correction
#reflectance.green <- ToTOAReflectance(landsat8=l, band="green", is.suncorrected = TRUE)
#Modified Normalized Difference Water Index (MNDWI)
#MNDWI <- ToMNDWI(l)
#The Normalized Difference Vegetation Index (NDVI), 
#Modified Normalized Difference Water Index (MNDWI) and Land Surface Water Index (LSWI) 
#need to compute the TOA Reflectance
#l <- ReadLandsat8(product,is.suncorrected = TRUE)
#MNDWI <- ToMNDWI(l)

（追記）

ToMNDWI関数をみたら

ToMNDWI
function (landsat8)
{
green <- ToTOAReflectance(landsat8, “green”)
swir1 <- ToTOAReflectance(landsat8, “swir1”)
mndwi <- (greeb - swir1)/(green + swir1)
return(mndwi)
}

greenがgreeb になっていた。

直っているコードがありました。
ToMNDWI.R