R

ランドサット7号のデータ2

raster ,cluster, som パッケージ

全然分からないけど土地被覆分類してみた

作業フォルダをtottori_city.tifのあるフォルダに変更

1
2
3
4
5
6
7
#setwd("tottori_city.tifのあるフォルダ")
library(raster)
library(RColorBrewer)
A_crop<- stack("tottori_city.tif")
#欠損値がないのを確認
summary(A_crop)  
DF <- as.data.frame(A_crop)

bandの情報を知る。（例）バンド２

1
2
3
4
5
6
7
8
9
10
A_crop[[2]]
# some basic statistics using cellStats()
cellStats(A_crop[[2]], stat=max)
cellStats(A_crop[[2]], stat=mean)
# This is equivalent to:
maxValue(A_crop[[2]])
# what is the maximum value of all three bands?
max(c(maxValue(A_crop[[1]]), maxValue(A_crop[[2]]), maxValue(A_crop[[3]])))
# summary() is useful function for a quick overview
summary(A_crop[[2]])

結果は省略

土地被覆分類(3パターン)

12種類に分類

k-means によるクラスタリングを行う
結果は、最初のクラスタのランダムな割り振りに大きく依存する。

1
2
3
4
5
6
7
E <- kmeans(DF, 12) #遅い
clusters <- raster(A_crop)   ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, E$cluster) 
clusters
#png("Classification01.png",width=1000,height=1000)
plot(clusters,col=brewer.pal(12, "Paired"))
#dev.off()

1
2
3
4
5
6
7
8
9
#cluster::claraを使う(kmeans関数を使うよりかなり速い)
library(cluster)
E=clara(DF,12)
clusters <- raster(A_crop)   ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, E$clustering) 
clusters
#png("Classification02.png",width=1000,height=1000)
plot(clusters,col=brewer.pal(12, "Paired"))
#dev.off()

自己組織化マップによるクラスタリング

1
2
3
4
5
6
7
8
9
library(som)
E.som<-som(DF, xdim=12, ydim=1)
#E.som$visual[,1]
clusters <- raster(A_crop)   ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, E.som$visual[,1]+1)  #番号を1から12にするために1加える
clusters
#png("Classification03.png",width=1000,height=1000)
plot(clusters,col=brewer.pal(12, "Paired"))
#dev.off()

1
2
3
4
5
#ゾーン統計（Zonal Statistics）
E.mean <- zonal(A_crop, clusters, fun="mean")  
E.min <- zonal(A_crop, clusters, fun="min")
E.max <- zonal(A_crop, clusters, fun="max")
E.sum <- zonal(A_crop, clusters, fun="sum")

結果は省略

放射輝度への変換方法を見つけたのでやってみた。

DNから放射輝度(Radiance)を求める方法Landsat-7 ETM+

デジタルカウント値を放射輝度に変換するには、各衛星データの仕様が必要。

p111r035_7x20001019.met

ダウンロードしたデータの必要な部分をコピペしてRに入力

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
#使用するデータによって異なる
LMAX_BAND1 = 191.600           
LMIN_BAND1 = -6.200          
LMAX_BAND2 = 196.500          
LMIN_BAND2 = -6.400           
LMAX_BAND3 = 152.900          
LMIN_BAND3 = -5.000    
LMAX_BAND4 = 157.400          
LMIN_BAND4 = -5.100          
LMAX_BAND5 = 31.060           
LMIN_BAND5 = -1.000         
#LMAX_BAND61 = 17.040           
#LMIN_BAND61 = 0.000           
#LMAX_BAND62 = 12.650           
#LMIN_BAND62 = 3.200           
LMAX_BAND7 = 10.800           
LMIN_BAND7 = -0.350         
#LMAX_BAND8 = 243.100         
#LMIN_BAND8 = -4.700          
##
QCALMAX_BAND1 = 255.0            
QCALMIN_BAND1 = 1.0             
QCALMAX_BAND2 = 255.0            
QCALMIN_BAND2 = 1.0               
QCALMAX_BAND3 = 255.0              
QCALMIN_BAND3 = 1.0            
QCALMAX_BAND4 = 255.0              
QCALMIN_BAND4 = 1.0           
QCALMAX_BAND5 = 255.0             
QCALMIN_BAND5 = 1.0           
#QCALMAX_BAND61 = 255.0             
#QCALMIN_BAND61 = 1.0            
#QCALMAX_BAND62 = 255.0              
#QCALMIN_BAND62 = 1.0            
QCALMAX_BAND7 = 255.0             
QCALMIN_BAND7 = 1.0           
#QCALMAX_BAND8 = 255.0             
#QCALMIN_BAND8 = 1.0

デジタルカウント値(DN)を放射輝度(L)に変換

計算式は
Radiance = ((Lmax-Lmin)/(QCALmax-Qcalmin)) * (QCAL-QCALmin) + Lmin

1
2
3
4
5
6
7
8
9
10
L1=((LMAX_BAND1 - LMIN_BAND1)/(QCALMAX_BAND1 - QCALMIN_BAND1))*(A_crop[[1]]-QCALMIN_BAND1)+ LMIN_BAND1
L2=((LMAX_BAND2 - LMIN_BAND2)/(QCALMAX_BAND2 - QCALMIN_BAND2))*(A_crop[[2]]-QCALMIN_BAND2)+ LMIN_BAND2
L3=((LMAX_BAND3 - LMIN_BAND3)/(QCALMAX_BAND3 - QCALMIN_BAND3))*(A_crop[[3]]-QCALMIN_BAND3)+ LMIN_BAND3
L4=((LMAX_BAND4 - LMIN_BAND4)/(QCALMAX_BAND4 - QCALMIN_BAND4))*(A_crop[[4]]-QCALMIN_BAND4)+ LMIN_BAND4
L5=((LMAX_BAND5 - LMIN_BAND5)/(QCALMAX_BAND5 - QCALMIN_BAND5))*(A_crop[[5]]-QCALMIN_BAND5)+ LMIN_BAND5
#L7の場合  A_crop[[6]]になるので注意 band=6は使っていないため
L7=((LMAX_BAND7 - LMIN_BAND7)/(QCALMAX_BAND7 - QCALMIN_BAND7))*(A_crop[[6]]-QCALMIN_BAND7)+ LMIN_BAND7
L <- stack(L1,L2,L3,L4,L5,L7)
A_crop
L

A_crop

class : RasterStack
dimensions : 302, 351, 106002, 6 (nrow, ncol, ncell, nlayers)
resolution : 28.5, 28.5 (x, y)
extent : 419990.2, 429993.8, 3926003, 3934610 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names : tottori_city.1, tottori_city.2, tottori_city.3, tottori_city.4, tottori_city.5, tottori_city.6
min values : 40, 18, 1, 1, 1, 1
max values : 255, 254, 255, 255, 255, 255

L

class : RasterStack
dimensions : 302, 351, 106002, 6 (nrow, ncol, ncell, nlayers)
resolution : 28.5, 28.5 (x, y)
extent : 419990.2, 429993.8, 3926003, 3934610 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=53 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names : tottori_city.1, tottori_city.2, tottori_city.3, tottori_city.4, tottori_city.5, tottori_city.6
min values : 24.170866, 7.179921, -5.000000, -5.100000, -1.000000, -0.350000
max values : 191.6000, 195.7012, 152.9000, 157.4000, 31.0600, 10.8000