Windowsがまだ3.1だったころ、某大学の学生時代、研究室でリモートセンシングの研究をしていました。衛星画像を解析して土地の被覆状況を調べるようなことをしていました。
これを読んでいる方の中にはまだ生まれていない方も多くいると思います。その当時はフロッピーディスクの時代で使えるPCも今では絶滅危惧種か絶滅しているようなものばっかりの中で、衛星画像の解析をしていました。最近Pythonを使って同じことを今の技術で行うと、いかに技術革新が進んだかわかったのでそのあたりをご紹介したいと思います。
平成6年
当時のIT周りの背景と使ったもの
この研究を始めたのは1994年、平成6年で元号が平成になって10年もたっていないころです。大ヒットしたWindows95がまた発売されていないころで、インターネットも一般では使われておらず大学などの研究で使われていたような時代です。ブラウザーは今はなきNetscape Navigatorの最初のバージョンがリリースされたころで、GAFAのうちAppleだけが存在していてAmazonは立ち上げたばかりのころです。また当時はフロッピーディスクの時代で、データのやり取りにCDすら使えなかった時代です。
当然、プログラミング言語についても今では使われないようなものばっかりで、Pythonはありましたが今のように機械学習で多く使われるようなものでもなく、そもそも機械学習という言葉自体あまり知れ渡っていないころです。自分も教師なしクラスター分析という言葉は研究中多く使ってきましたが、機械学習という言葉を使ったことはあまり記憶にありません。
そのような環境で当時大学の研究室で使えたものが以下の通りです
- Windows 3.1 32bit
- 5インチ、3.5インチフロッピーディスク（容量は3.5インチで1.4Mほどです）
- RAM 4-64M, 1cpu
- Visual Basic2.0
今IT業界で仕事をしていますが、これはほぼ化石レベルの古さですね。
衛星画像データ
使用した衛星画像データはランドサット5号のTMデータで、研究費用で日本各地のデータを購入しました。前述のように一般に使われている保存媒体は1.4M のフロッピーディスクでインターネットもまともに利用されていないころですので、当然衛星画像データをダウンロードなど出来ません。
購入したデータはフィルム映画の映写機にのような機械を使用して読み取る磁気テープで、何時間もかけてテープからデータを読み取り、保存したデータを使用するパソコンへコピーといった作業をしていました。
使用したデータは以下の通りです
- ランドサット5号: 1992年4月21日
- 分解能: 30m （30メートル以上のものを見分けることとができます。）
- 解析都市: 京都市近郊
- 解析範囲: 縦501、横325ピクセル（山崎ジャンクションから左京区北白川ぐらいの範囲を長方形で切り取ったぐらいです）
- 使用バンド: 1-7
教師なしクラスター分析
当時の環境では、ほぼ選択肢がなかったように記憶していますが、非階層クラスター分析を使っています。教師あり、なしについて同じ研究室で教師ありを行っている人がいたので教師なしを使うことにしたはずです。また当時は今のようにGoogle Mapなどなかったので、教師ありにすると土地被覆状況をあるていど実測しないといけなくなるので、教師なしにしたように思います。基本的なロジックはISODATA (Interactive Self-Organization of Data)法を少し変更したものをもとにVisual Basic2.0を使って開発しました。
Visual Basic2.0で機械学習
Pythonのようにライブラリーがなかったので、基本的には自力でアルゴリズムを作る必要があります。
パラメーターの決定
初期クラスタの数、再配置の収束条件、解析の終了条件などは解析を始めるまえに入力してあらかじめ設定しておきます。
| 項目 | 値 |
|---|---|
| 使用バンド | BAND1-7 |
| 初期クラスタ数 | 7 |
| 再配置収束条件 | 5% |
| 解析終了条件 | 0.7% |
初期クラスタの重心の決定
初期クラスタ数: n、バンド: b、クラスタ: c
として各バンドの平均値αと標準偏差βをもとに適当な重心値を計算しています。
^{b} g _c = ^{b} \alpha _c + ^{b} \beta _c
\left( \frac{ 2\left(c -1 \right)}{n-1} - 1\right)
For c = 1 To n
g(b,c) = a(b) + std(b) * ((2*(c-1))/(n-1) -1)
Next c
再配置法
バンド: b、クラスタ: c、使用バンド数： k
それぞれの地点μと各クラスタの重心χとのユーグリッドの距離を計算し、最も近いクラスタに再配置をします。
d _c =
\left\{ \sum_{b=1}^{k} \left(^{b} \chi_c - ^{b} \mu \right)^2\right\}^\frac{1}{2}
'各個体とクラスターの重心との距曜の計算
For yl = starty To endy
For xl = startx To endx
min = 10^4 '初期最低距離
place = XMAX*(y1-1) + x1 + head
For c = 1 To n
For b = 0 To k-1
Get filenum(b), place,resdata
rdata = Asc(resdata)
d1 = (g(b,c)- rdata)^2
d(c)= d (c) + d l
Next b
'距離が最低距離より小さい場合は再配置して最低距離の再設定
If d(c)< min Then
min = d(c)
classnum = c
End If
d(c) = 0
Next c
Next xl
Next yl
最大64Mのメモリーではデータを一括で読み取ってメモリーにキープできるわけもないので、ループでいちいち読み取っています。
再配置収束条件
再配置後に所属するクラスタを変えた個体数が再配置収束条件5%以上の場合は、再配置後のクラスタの重心を再計算し再配置を繰り返します。5%以下の場合は収束したとみなします。
'クラスタを変えた個体数算出
Get filenumB,place,Gcmoji
Gclassnum = Asc(Gcmoji)
If Gclassnum <> classnum Then
change = change + 1
'クラスタ番号を分析結果ファイルに書き込み
Pcmoji = Chr$(classnum)
Put filenumB,place,Pcmoji
End If
'収束条件
calend = 5
Pchange = 100 * (change / X*Y)
If Pchenge <= calend Then
MsgBox "計算が終了しました"
End If
解析終了条件
再配置収束後、クラスタの各バンドごとの標準偏差を求め、最大値を持つクラスタを割り出します。最大値が、画像全体のそのバンドの標準偏差と比較して70%以上の時はクラスタを2分割します。標準偏差が大きすぎるクラスタを減らす作業をしています。
標準偏差が大きすぎるクラスタはこんな感じです。
これをそのままにしておくとこのクラスタの土地被覆状況を判断するのは難しくなります。
'標準偏差が最大値のクラスタをバンドごとに求める
For b = O To k-1
stdmax(b) = 0
For c = l To n
std(b, c)= Sqr(gg(b, c)-g(b, c)^2)
If stdmox(b) < std(b, c) Then
stdmax(b) = std(b, c)
stdmac(b) = c
End If
Next c
Next b
'標準偏差が最大のクラスタを分裂し、クラスタ番号を追加して重心を再計算
For b = O To k-1
If stdmax(b) > stdall(b) * 0.7 Then
np = np + 1
For fc = stdmac(b) To n + np Step n + np - stdmac(b)
For fb = O To k-1
g(fb,fc) = GetCenter(fb,stdmac(b))
Next fb
Next fc
End If
Next b
'クラスタ総数の追加
n = n + np
解析結果
これら以外にも、各クラスタに色を割り当てて画像を保存するコードや、各クラスタの個体数を計算するプログラムなどほぼほぼすべて手書きでロジックを考えて行った結果が、以下の通りです。
当時の研究のコピーで白黒しかないのであまりはっきりとは出ませんが、桂川、鴨川、宇治川などがまずまずはっきりと出ています。また御所や二条城のあたりも、それがあるというのがわかります。東山のところは山と街の境がはっきりと出ていて、土地の被覆状況が異なるのが分かります。
| 項目 | 値 |
|---|---|
| 計算回数 | 22 |
| 計算開始日 | 12月26日20:52:33 |
| 計算終了日 | 12月28日14:35:17 |
| 分裂したクラスタ数 | 8 |
見ての通り22回の反復計算をするのに40時間以上かかっています。
分裂したクラスタ数が８なので、合計クラスタ数は15となります。しかしながら、実際クラスタを分類すると特に小さいクラスタはどのように土地が使われているかはっきり分らず、未分類となってしまったクラスタが4つ出てきました。また、同じ用途で使われていても、違うクラスタとして判別されたものもいくつか出てきました。現在のようにGoogle Mapなどがあればもう少し正確な分類ができたかもしれませんが、当時は地図と比べての分類だったので、これぐらいが限界でした。その内訳は以下の通りです。
| クラスタ番号 | 個体数 | 分類 |
|---|---|---|
| 1 | 2017 | 水域 |
| 2 | 9040 | 森林 |
| 3 | 3045 | 未分類 |
| 4 | 27303 | 住宅地 |
| 5 | 14798 | 水田 |
| 6 | 14836 | 住宅地 |
| 7 | 1618 | 未分類 |
| 8 | 5404 | 未分類 |
| 9 | 7103 | 森林 |
| 10 | 7374 | 荒地 |
| 11 | 27081 | 市街地 |
| 12 | 4209 | 草地 |
| 13 | 16946 | 市街地 |
| 14 | 13406 | 水田・畑 |
| 15 | 9472 | 未分類 |
その結果、解析範囲の京都市近郊の土地被覆状況は以下の通りでした。
| 土地被覆分類 | 割合 |
|---|---|
| 住宅地 | 25.7% |
| 市街地 | 24.7% |
| 水田・畑 | 10.3% |
| 森林 | 9.8% |
| 水田 | 9.0% |
| 荒地 | 4.5% |
| 草地 | 2.5% |
| 水域 | 1.2% |
| 未分類 | 11.9% |
令和3年
解析に使用したもの
最新ではありませんが当然平成6年のころに比べたら格段に良くなっています。
- Windows10
- 64bit
- RAM 16G, 2 cores, 4プロセッサ
- Python3.8
衛星画像データのダウンロード
20数年前とちがい現在はインターネットの環境も整っていますし、利用規約範囲内であれば無料でダウンロードして使用することもできます。
【実習編】～Landsat8衛星画像で植生を映えさせた地図を描こう～ を参考に、https://landbrowser.airc.aist.go.jp/landbrowser/index.html より以下のデータをダウンロードして、今回の解析に使いました。
- ランドサット8号: 2020年10月27日
- 分解能: 30m （30メートル以上のものを見分けることとができます。）
- 解析都市: 京都市近郊
- 解析範囲: 縦500、横300ピクセル（山崎ジャンクションから左京区北白川ぐらいの範囲を長方形で切り取ったぐらいです）平成6年と完全に一致はしていませんが、ほぼ同じ個所を同じようなサイズで切り取りました。
- 使用バンド: 1-7
k-meansとx-meansで解析
パラメーターの決定
使用するバンド数は前回同様7です。ライブラリーを使用するので、計算終了条件などは指定する必要がありません。
ただ、k-meansの場合はクラスタ数をあらかじめ決めている必要があり、x-meansの場合は最大クラスタ数を設定しないとクラスタ数が多くなりすぎ、最終的に土地の被覆状況が分類できなる恐れがあります。
まず、エルボー法でクラスタ数を求めてみました。
def k_get_elbow(X,cluster_num):
distortions = []
for i in range(1,cluster_num+1):
km = cluster.KMeans(n_clusters=i)
km.fit(X)
distortions.append(km.inertia_)
#y_km = km.fit_predict(X)
plt.plot(range(1,cluster_num+1),distortions,marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.show()
結構なだらかな曲線で表示され、はっきりとどのあたりが一番良いかというのは、今回使用したデータではわかりませんでした。
また、x-meansで最大クラスタ数を例えば30ぐらに指定すると、30個までクラスタ数が分けられるような結果になります。クラスター分析という意味ではおそらく良い結果なのでしょうが、今回の目的である土地被覆状況を分類し以前の結果と比較する目的ではかなり使いにくい結果となりました。コードは後程紹介しますが、最大クラスタ数を30にしてx-meansで解析した結果はこのような感じです。
ほぼ、地図みたいになって、詳細すぎてどのクラスタが何を表しているかを判別するのはほぼ不可能です。
したがって、今回は前回の解析結果と比較するのが目的なので、解析の精度という点ではあまり深く追求せずk-meansのクラスタ数、x-meansの最大クラスタ数ともに10として解析を進めました。
| 項目 | 値 |
|---|---|
| 使用バンド | BAND1-7 |
| k-means クラスタ数 | 10 |
| x-means 最大クラスタ数 | 10 |
データの読み取り
以前と違ってメモリーも十分あるのでデータを一気に読み取ります。データは各バンドごとに読み取り7次元のndarrayを返します。
def k_get_raw_data(city,dataset,r,x,y,h,w):
mimage = np.zeros((h, w, r), dtype=int)
for i in range(r):
b = i + 1
with rasterio.open(city + '\\' + dataset + str(b) + '.tiff') as src:
data = src.read()
img = data[0][y:y+h, x:x+w]
mimage[:, :, i] = img
new_shape = (mimage.shape[0] * mimage.shape[1], mimage.shape[2])
X = mimage[:, :, :r].reshape(new_shape)
return X
k-meansで解析
前述のVisual Basicの場合はデータの読み取り、再計算、分裂など多くの行数を費やしてロジックを作りました。紹介したものもほんの一部で、もっと長いコードを使っています。（自分のプログラミング技術も相当乏しかったのですが）
k-meansを使って解析すると解析部分は数行で済ませることができます。先ほど読み込んだデータと、初期設定のクラスタ数をcluster.KMeans(n_clusters=cluster_num)に渡すだけです。
from sklearn import cluster
def k_cal_show_image(X,cluster_num):
k_means = cluster.KMeans(n_clusters=cluster_num)
k_means.fit(X)
解析結果は最終的には先ほどのように画像に表示したいので、それぞれの地点のクラスタ番号を読み取って新たに画像サイズにあったndarrayを作ります。
X_cluster = k_means.labels_
X_cluster_new = X_cluster.reshape([h,w])
x-meansで解析
k-meansと同様に解析部分は数行です。
from pyclustering.cluster.xmeans import xmeans, kmeans_plusplus_initializer
def x_cal_plot(X,cluster_num):
# 初期クラスタ数を２からスタート
initial_centers = kmeans_plusplus_initializer(X, 2).initialize()
analysis = xmeans(X, initial_centers, kmax=cluster_num, tolerance=0.025, criterion=0, ccore=True)
analysis.process()
こちらも解析気結果を画像で表示するために、各クラスタにクラスタ番号を割り当てそれぞれの地点に書き込みします。
例えばあるクラスタにクラスタ番号７を割り当てて、それに属する個体の位置を表示すると以下のようになります。
print('Cluster ' + str(7) + ' first 10: ' + str(clusters[7][:10]))
Cluster 7 first 10: [593, 605, 769, 1128, 3284, 3285, 3556, 3581, 3582, 3769]
これは593番目、605番目・・・はクラスタ７に所属するということです。
クラスタ番号の割り当てを行って、先ほど同様に画像サイズにあったndarrayを作ります。
clusters = analysis.get_clusters()
X_cluster_new = np.zeros(X.shape[0], dtype=int)
for i in range(len(clusters)):
for j in range(len(clusters[i])):
X_cluster_new[clusters[i][j]] = i
X_cluster_new = X_cluster_new.reshape([h,w])
画像保存
k-means、x-meansともに各クラスタに任意の色を割り当てて解析結果を画像として保存します。
import matplotlib.patches as mpatches
im = plt.imshow(X_cluster_new,cmap='jet')
colors = [ im.cmap(im.norm(value)) for value in cluster_size[:,0]]
patches = [ mpatches.Patch(color=colors[i], label="Cluster {l}".format(l=cluster_size[:,0][i]) ) for i in range(len(cluster_size[:,0])) ]
plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0. )
各クラスタの割合を表す円グラフも表示します。
plt.title('% of clusters')
plt.pie(x=cluster_size[:,1], labels=cluster_size[:,0], autopct='%.2f%%',colors=colors)
解析結果
k-meansとx-meansを使った解析結果を比較します。
k-means
x-means
色合いは異なりますが、出された結果はかなり似通っていました。平成6年の結果と違いカラーで表示されているのでさらにわかりやすくなっています。桂川、鴨川、宇治川はもちろん、京都競馬場、二条城、御所などもはっきりと出ています。また、名神高速道路や、平成6年にはなかった第二京阪道路なども見て取れます。二条城や御所には多くの木々があるのですが、それが東山の森林などと同じクラスタに分類されているのもわかります。また、街中の碁盤の目の道路も大通りの分は読み取ることができます。 （地名などを入れた画像は後の解析画像の比較をご覧ください）
変動係数(Coefficient of variation)
変動係数(Coefficient of variation)を比較してみましたが両方とも下のグラフのようにほぼ同じ感じで出ました。
これはx-meansでの変動係数ですが、k-meansの結果もほとんど同じでした。同じクラスタでも、Band5で係数が大きくなっているのが分かります。これは、次のクラスタごとのヒストグラムでも見て取ることができます。
ヒストグラム
ヒストグラムを作ってばらつきの少なかったクラスターと大きかったものを比べてみました。
まずは、先ほどの変動係数のグラフで係数が小さかったクラスタ３と４です。全体的にきれいな山型になっていますが、変動係数のグラフ通りBand5でばらつきが大きくなっています。
それに比べてばらつきが最も大きかったクラスタ6と7です。クラスタ7はばらつきが大きすぎるので変動係数のグラフからは省いています。全体的に偏っていたり、山型がきれいに出ていませんし飛び地もあったりします。
ばらつきが目立ったのがBand5で、これは植生を判別するのにつかわれます。下のクラスターの分類を見ての通り、多くのクラスターで植生が入っているのでこれが原因かと考えています。また、クラスタ6と7は小さいクラスタで、どうしてもばらつきが出やすかったと考えられます。
結果
| 項目 | 値 |
|---|---|
| 計算回数 | 55(k-meansの場合) |
| 計算時間 | k-means、x-meansともに1分 |
画像や変数係数を比較してもk-meansとx-meansでほぼ差がなかったので、土地被覆分類は先ほど変数係数を用いたx-meansの結果を元に分類しました。
| クラスタ番号 | 土地被覆分類 | 割合 |
|---|---|---|
| 3 | 住宅地 | 28.60％ |
| 4 | 市街地 | 28.47％ |
| 2 | 畑・水田 | 8.96％ |
| 1 | 芝、野原 | 8.79％ |
| 8 | 森林 | 8.34％ |
| 0 | 水田・畑 | 5.91％ |
| 9 | 芝・広葉樹 | 5.34％ |
| 5 | 水/水辺/芝 | 4.61％ |
| 6 | 工場・ビル | 0.91％ |
| 7 | 工場・ビル | 0.06％ |
平成6年と令和3年の解析結果の比較
平成6年にWindows3.1とVisual Basic2.0でほぼ自前でロジックを組んで解析した結果と、令和3年にWindows10とPython3.8でライブラリーを活用して解析した結果を比較します。比較・考察するのは解析のIT技術的なところと、実際の土地被覆状況などの2点です。
IT技術的な比較
- 一番の違いはやはり計算時間の差です。平成6年の解析では40時間以上かかったのが令和3年の解析は1分で終了しています。
- 平成6年の解析と令和3年のk-meansの解析は同じユーグリッドの距離を使っているのですが、平成6年は22回の計算を40時間以上かけているのに対して、令和3年は55回の計算を1分で終わっています。
- 当然ですが、使っている環境が大きく違いメモリーは最大64Mと16G、プロセッサも数だけでも1と4の違いもありますし、性能もはるかに良くなっています。データのI/Oではいちいち読み書きしていたのが、すべてメモリー内で処理され最終結果のみ書き出しているのでI/Oに費やされる時間も大きな差が出ているはずです。
- プログラミングでは、Visual Basic2.0でロジックをすべて手書きしているのでかなりの行数を使います。ソースコードをすべて保管していないのですが何ページ分もあったのを覚えています。一方、Pythonだとライブラリーを使うので、全体的にはるかに行数も少なくすんでいます。
- 画像処理以外のデータの処理（変動係数、ヒストグラム、クラスタ数など）も、使ったライブラリーをもとに変数を取り出すだけなど簡単な処理で、新たにロジックなどを作ることもなくできました。
- ライブラリーで簡単にできてしまうので、統計学的な知識がなくてもある程度のプログラミングができれば誰でも簡単にできるのではと思いました。AIや機械学習を軽く勉強するには、かなり取り掛かりやす物だと思います。また、平成6年の解析ではクラスタの精度を高めるために標準偏差を用いて分裂などの作業をしていましたが、令和3年の解析ではライブラリーですべて行ってくれてヒストグラムでも見て取れる通りある程度の精度のクラスタを手間をかけることなく作れました。
- 解析用のデータを取得するのも、平成6年はデータをテープ状態で購入してそこから読み取るといた作業をしていましたが、今回はサイトよりダウンロードしてすぐに使用することができました。一見普通のことのようですが、二十数年前では考えられなかったことです。
土地被覆状況の比較
解析画像の比較
考察
- 目立った違いは平成6年の解析では未分類が約12%も出たことです。プログラムの精度、ランドサットの違いもありますが、大きな原因の一つは判別方法だと思っています。当時はGoogle Mapなどないですし現地に向かう時間もなかったので、地図などをもとに分類を判別しました。紙の地図とGoogle Mapの衛星写真をつかって判別するのでは大きく精度も異なると思っています。未分類になったのは、地図と比較してはっきり何に使われているかわからなかったところで、市街地や住宅地のようにはっきりとわかるところではなく、水田と畑、芝地と荒地、森林と草地など分かりにくかったところが大半だっと記憶しています。
- 上の二つの解析画像を比べるとわかるのですが、1992年の画像は2020年の画像に比べて若干反時計回りに傾いています。今回はそれを是正するような詳細な解析はしておりませんが、これにより2020年の解析範囲のほうが東山近辺の森林地区が多く含まれ、向日市・長岡京市の市街地範囲が少なく含まれています。
- 平成6年と比べて令和3年解析では市街地・住宅地が約7%ほど増加しています。使用したランドサットは共に30メートルの分解能なので約9.4㎢ほど市街地・住宅地が1992年から2020年の28年ほどの間で増えたということでしょうか。それとは逆に水田や畑は5%ほど減少しています。前述の通り、平成6年にはなかった第二京阪道路なども見て取れます。これは市街地として分類されたクラスタ4に含まれており、傾きによって含まれる市街地部が少ないにもかかわらず市街地が増加したのは、このような都市開発が市街地・住宅地の増加に関連していると考えられます。
- 平成6年の解析では森林は1つのクラスタで分類していますが、今回は森林部分と広葉樹や芝を含んだ部分の2つのクラスタで分類できています。クラスタ９が広葉樹、芝部分なのですが桂川・宇治川河岸の樹木のあるところがはっきりと出ています。ただ、近辺にある芝部分との区別はできていないようです。
- 令和3年の解析では、前回なかった工場・ビルと分類されたクラスタが出てきました。Google Mapと見比べてもはっきりとそれとわかるクラスタ6及び7が、山崎ジャンクションおよび太秦の近辺に見られます。京都駅もクラスタ6で分類できています。
- 平成6年の解析では御所や二条城があるというのはわかる程度でしたが、令和3年の解析ではそれぞれの樹木の部分とそうでないところがはっきりと分類できています。
- 両方とも宇治川の南の巨椋近辺の水田、畑はかなり明確に表示されています。ある程度大きな範囲で被覆状況が同じであれば、27年前のプログラムでもそれなりに分析を行えていたのが分かりました。当時の自分の技術力を考えると、悪くないのではと思っています。
京都競馬場拡大図
- クラスタではありませんが、令和3年の解析では京都競馬場もはっきりと見分けることができます。外の芝のコースは芝・広葉樹と分類されたクラスタ9、ダートのコースは畑・水田と分類されたクラスタ2となっています。10月末のデータなので実際の畑・水田は水も抜けて土状になっていると思われ、それとダートコースが同じクラスタと分類されたのではと考えられます。また、京都競馬場の中心部は池なのですが、水/水辺/芝と分類されたクラスタ5と分類されています。
その他の解析
拡大範囲で解析
次に、比較をするために解析範囲を少し拡大して分析をしてみました。使用しているのは全く同じランドサットデータと解析プログラムです。
| クラスタ番号 | 土地被覆分類 | 割合 |
|---|---|---|
| 6 | 市街地 | 25.62％ |
| 7 | 住宅地 | 17.77％ |
| 1 | 森林・西向き | 12.67％ |
| 0 | 森林・東向き | 10.07％ |
| 2 | 広葉樹 | 9.70％ |
| 5 | 水田 | 9.37％ |
| 8 | 畑 | 6.18％ |
| 3 | 水/水辺 | 4.10％ |
| 4 | 芝、野原 | 3.89％ |
| 9 | 工場・ビル | 0.66％ |
全体の数が増えたせいか、機械学習がより正確にクラスタを分類できるように判別できたと思います。森林の中での違いやと広葉樹林の違いなどを判別できたり、河岸も河川と岸をはっきり区別できています。また、市街地でも高速道路などがよりはっきりと表れました。
京都は盆地と言われますが、この画像からもはっきり市街地住宅分が東西および北側の山に囲まれた盆地であることが良くわかります。
温度分布
京都は盆地で夏は暑く冬は底冷えと言われますが、実際どんな感じか、同じ市内で差があるのか、それぞれの時期のランドサットデータで比べてみました。熱赤外はバンド 10、11で使用され、それを使って解析すると地表面温度を見ることができます。それを使って2017年の12月と翌年7月の比較をしてみました。
赤色が濃いのは温度が高く、薄くなるほど温度が低いことを表しています。
緑で囲まれた御所の西の部分7月は特に色が濃くなっています。同じ場所を12月の左側の画像と比較すると、他の地域より薄くなっています。それと比較して黄色で囲まれた右京区の嵐山から西院近郊は夏は、御所西側より色が薄く、冬には濃くなっています。
このことからすると、御所西側は他の地区と比較して夏は暑く冬は寒いと考えられます。逆に嵐山から西院近郊は、同じ京都市内でも比較的夏は涼しく冬は暖かく住みやすいのではとも考えられます。2日分のデータなので断言できるようなものではないですが、機械学習をつかうとこのようなことも調べることができるとわかりました。
まとめ
- 27年間でITの世界は全く違うものになりました。この業界がいかに早く進化しているか、改めて考えられるきっかけになりました。
- 27年前の技術でもロジックさえしっかりしていれば、ある程度の解析ができていたことが判明しました。当然今の技術を使った方が、精度もよく時間も短くできますが、あの当時のレベルでもそれなりの解析ができていたのは良かったと思っています。
- 20年以上前のプログラムは今はほぼ役に立つことはないですが、クラスター分析をするための統計学的な基礎は一般的に機械学習と言われるITの技術の中でも生かせることが分かりました。
- 今回のような簡単な解析でも、リモートセンシングを使って土地の被覆状況の変化を見て取ることができました。市街地・住宅地が広がっていることや、水田・畑が減少していることなど正確でないにしてもより精度を高くすれば環境問題などにも利用できるものと思われます。
- 温度分布を解析することによって、その地域で比較的住みやすい気候の所とそうでないところを見つけることができそうです。
参考文献および注釈
他にも多くのサイトを参考にさせていただきましたが、主たるものは以下の通りです。
- 画像解析ハンドブック 単行本 – 1991/1/1 高木幹雄 (著), 下田陽久 (著)
- リモートセンシング
- Netscape Navigator
- フロッピーディスク
- ランドサット5号
- ランドサット8号
- 【実習編】～Landsat8衛星画像で植生を映えさせた地図を描こう～
- ランドサットデータについて (https://landbrowser.airc.aist.go.jp/landbrowser/index.html)
- “ The source data were downloaded from AIST’s LandBrowser, (https://landbrowser.airc.aist.go.jp/landbrowser/). Landsat 7/8 data courtesy of the U.S. Geological Survey.”