Go to Qiita Advent Calendar 2024 Top
LoginSignup
12
15

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@dg4101

熱海土石流災害を衛星画像と機械学習で解析

Last updated at Posted at 2021-08-26

2021年7月に起こった熱海での土石流災害を衛星画像と機械学習を用いて分析してみました。

以前、下の記事で京都の土地被覆状況を衛星画像、機械学習とクラスター分析を用いて解析する方法を書いてみましたが、今回はそれを応用してより、実際起こっている環境問題などに取り組んでみました。

#衛星画像データのダウンロード
今回の解析では、Sentinel2のデータを使って分析をしてみました。ダウンロードする方法などは、参考文献および注釈にその方法を説明している記事をリンクしましたのでそちらを参考にしてください。

  • Sentinel2 :　
    • 2021年8月5日　土石気流発生後、晴天で雲が少ない日のデータ
    • 2020年8月15日　土石気流発生前と比較するため1年前の同時期の晴天で雲が少ない日のデータ
  • 分解能: 10m　（10メートル以上のものを見分けることとができます。）
  • 解析都市:　熱海近郊
  • 使用バンド:　2、3、4、8。10mの分解能のバンドのみを使用しています。

#教師なしクラスター分析

##解析範囲の切り出し
方法はいろいろとあると思いますが、今回は画像を表示してそこで解析範囲を指定しました。Sentinel2のデータにはTCI(True Color Image)が含まれていますので、それを使うとはっきりとした画像で解析範囲を探すことができます。

import numpy as np
import matplotlib
import rasterio

# 画像データ
with rasterio.open('T54SUD_20210805T012701_TCI.jp2') as src2:
    data2 = src2.read()

plt.imshow(data2[0])
plt.show()

ズームボタンで必要な大きさまで拡大して、右下に表示される座標を元に切り出し範囲を指定します。
image.png
また黄色で囲んだ部分の明るい部分は雲で、これにより解析を行う熱海市街の海岸近くは雲がかかっていないことが見て取れます。

##x-meansで解析
x-meansを使って教師なしクラスター分析を行いました。x-meansの詳細に関しては参考文献および注釈から参考にさせていただいた記事をご覧ください。

###データの読み取り
前項の画像を元にまずは、災害地区を含んだ全体像の解析をします。
範囲としては縦1000ピクセル、横500ピクセルです。1ピクセル10mなので、10Kmx5Kmの範囲となります。
コードは京都の記事とほぼ同じですが、使用バンドを指定できるように改良しました。

'''
# 10m
x=2200
y=1000
h=1000
w=500
r =[ '02', '03', '04','08']
cluster_num = 15
'''
import numpy as np
import rasteri
def k_get_raw_data(city,dataset,r,x,y,h,w):
    mimage = np.zeros((h, w, len(r)), dtype=int)
    for i in range(len(r)):
        b = r[i]
        with rasterio.open(city + '\\' + dataset + b + '.jp2') as src:
            data = src.read()

        # 画像サイズの確認
        print('Shape: ' + str(data.shape))
        print('Dimension: ' + str(data.ndim))    
        print('Dataset:' + dataset + b )
        org_img = data[0]

        img = org_img[y:y+h, x:x+w]
        print('Extracted Image Shape: ' + str(img.shape))
        mimage[:, :, i] = img

    new_shape = (mimage.shape[0] * mimage.shape[1], mimage.shape[2])
    X = mimage[:, :, :len(r)].reshape(new_shape)
    return X

###機械学習
解析結果の保存部分でも、前回より少し改良してみました。

clusters.sort(key=len, reverse=True)

画像で表示した時に、クラスター番号順に色を表示できるようにしています。

np.save(imagefile + '.data',X_cluster_new)

座標ごとのクラスター番号の割り当てをndarrayデータとして後の解析のために保存しています。

import numpy as np
import matplotlib.pyplot as plt
from pyclustering.cluster.xmeans import xmeans, kmeans_plusplus_initializer

def x_cal_plot(X,city,r,cluster_num,h,w,now):
    # 画像ファイル
    listToStr = ' '.join(map(str, r))
    imagefile = 'images/x-' + city + '_bandnum_' + listToStr + '_clusternum_' + str(cluster_num) + '_' + now

    # 初期クラスター数２で開始
    initial_centers = kmeans_plusplus_initializer(X, 2).initialize()
    instances = xmeans(X, initial_centers, kmax=cluster_num, tolerance=0.025, criterion=0, ccore=True)
    instances.process()
    clusters = instances.get_clusters()
    clusters.sort(key=len, reverse=True)

    X_cluster_new = np.zeros(X.shape[0], dtype=int)

    for i in range(len(clusters)):
        print('Cluster ' + str(i) + ': ' + str(len(clusters[i])))

        X_cluster_info[0:,i] = i
        X_cluster_info[1:,i] = len(clusters[i])

        # Insert cluster num into X_cluster_new
        for j in range(len(clusters[i])):
            X_cluster_new[clusters[i][j]] = i

    X_cluster_new = X_cluster_new.reshape([h,w])

    #############################################
    # 画像
    fig = plt.figure(dpi=400)
    fig.suptitle('X-Means ' + city + ': cluster:' + str(len(clusters)) )
    plt.get_current_fig_manager().full_screen_toggle()
    plt.rcParams["font.size"] = 6

    import matplotlib.patches as mpatches
    im = plt.imshow(X_cluster_new,cmap='jet') 
    colors = [ im.cmap(im.norm(value)) for value in range(len(clusters))]
    patches = [ mpatches.Patch(color=colors[i], label="Cluster {l}".format(l=i)) for i in range(len(clusters)) ]
    plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0. )    
    plt.savefig(imagefile + '.png')
    np.save(imagefile + '.data',X_cluster_new)

#解析結果

##熱海全域

2020年8月15日 2021年8月5日
image.png image.png
2020年と2021年の解析結果を比較すると、赤で囲まれた部分に明らかな違いが出ています。熱海市 大規模土石流 現地調査 写真レポートを参考にすると、その範囲で土石流災害が起こっているのが分かります。また、海岸線に2020年ではなかった海水とは違うクラスターが見られます。これは、雨によって流れ出た土砂が原因で海水が濁っているからではと考えています。

##True Color画像の表示
Sentinel2のバンド2（青）、バンド3（緑）、バンド4（赤）を組み合わせてTrue Colorの画像を作って先ほどの土石流災害の箇所や海岸線がどの様になっているか見てみます。

import rasterio
from rasterio import plot
from osgeo import gdal

def get_color_pic(city,dataset,r,x,y,h,w,now):

    b = r[0]
    with rasterio.open(city + '\\' + dataset + b + '.jp2',driver='JP2OpenJPEG') as src0:
        data0 = src0.read()
        img0 = data0[0][y:y+h, x:x+w]
    b = r[1]
    print('image:' + str(b))
    with rasterio.open(city + '\\' + dataset + b + '.jp2',driver='JP2OpenJPEG') as src1:
        data1 = src1.read()
        img1 = data1[0][y:y+h, x:x+w]
    b = r[2]
    print('image:' + str(b))
    with rasterio.open(city + '\\' + dataset + b + '.jp2',driver='JP2OpenJPEG') as src2:
        data2 = src2.read()
        img2 = data2[0][y:y+h, x:x+w]

    # Set image file name
    listToStr = ' '.join(map(str, r))
    imagefile = 'images/x-' + city + '_bandnum_' + listToStr + '_' + now + '_color'


    trueColor = rasterio.open(imagefile + '.tiff','w',
                    driver='Gtiff', 
                    width=w, height=h, 
                    count=3, crs=src2.crs,
                    transform=src2.transform, 
                    dtype=src2.dtypes[0]) 
    trueColor .write(img0,3) # blue
    trueColor .write(img1,2) # green
    trueColor .write(img2,1) # red
    trueColor .close()

    scale = '-scale 0 255 0 25'
    options_list = [
        '-ot Byte',
        '-of JPEG',
        scale
    ] 
    options_string = " ".join(options_list)

    gdal.Translate(imagefile +'.jpg',
                imagefile +'.tiff',
                options=options_string)

    src = rasterio.open(imagefile + '.jpg')
    plot.show(src)

image.png
土石流が起こったところがかなり明確に表示されていますし、また海岸線も濁ったような色が出ています。

##土石流発生箇所の解析
次に土石流が発生した箇所に絞って解析をしました。コードは同じですが、解析範囲を縦160ピクセル、横200ピクセルです。1ピクセル10mなので、1.6Kmx2Kmの範囲となります。

'''
# 10m
x=2400
y=1165
h=160
w=200
'''
2020年8月15日 2021年8月5日 　クラスター
image.png image.png image.png
image.png image.png

1年前の土地被覆と比較しても、土石流が山から発生して町の中に流れ込んでいる様子がはっきりと表れています。True Color表示で見て取れるように土石流は町中まで流れて、解析画像から見ても土地被覆が変わっているのが分かります。

##土石流発生が流れた範囲の抽出
次に、土石流が流れた範囲がどれ位だったか調べてみました。方法としては先ほどの2つの解析画像で明らかにクラスター番号が大きく変化した箇所、つまり土地の被覆状況が大きく変わった個所を調べてみました。

import numpy as np
import matplotlib.pyplot as plt
from rasterio import plot
import rasterio
import cv2

# 読み込み
data_before = np.load("images/x-atami_bandnum_02 03 04 08_clusternum_15_20210814-1903_20200815_.data.npy")
data_after = np.load("images/x-atami_bandnum_02 03 04 08_clusternum_15_20210814-1903_20210805_.data.npy")

# クラスター番号の置き換え
data_load = np.uint8()
data_load = np.where((5 < data_after) & (data_after < 10) & (-1 < data_before) & ( data_before < 5), 20,data_after)
data_load = np.where(( data_after == 9) & (4 < data_before) & ( data_before < 9) , 20,data_load)
# 土石流個所を白色
data_load = np.where(( data_load != 20), (data_load + 1 ) * 15, 255)

# 白色のピクセル数
print(data_load.shape)
print("White pixel: " + str(np.count_nonzero(data_load == 255)))

# 画像の保存
data_load8 = (data_load).astype(np.uint8)
data_load_cv = cv2.applyColorMap(data_load8, cv2.COLORMAP_BONE)
cv2.imshow('image with colormap', data_load_cv)
cv2.waitKey()
cv2.imwrite('images/x-atami_bandnum_20210805_change_extract.png',data_load_cv)

2020年と2021年を比較し土石流で土地被覆が明らかに変わった個所のクラスター番号を変え、画像で保存するときに白色になるようにします。

  • 2020年でクラスター0-4の森林地区と分類されているが、2021年ではクラスター6-9となっているところ。
  • 2020年でクラスター5-9の住宅地区近郊と分類されているが、2021年ではクラスター9（おそらく土と分類されている）となっているところ。

置換が終わって保存した画像データは、以下の通りです。

image.png

土石流が発生したところが、白入りに代わっています。そのほかの箇所も変わっていますが、土石流がどれくらいの範囲で起こったか調べるため、今回は赤で囲ったところのみを調べます。

画像ツールで、赤で囲ったところ以外を切り取ります。
image.png
これで、土石流発生個所のみが白色で、それ以外は違う色で表すことができました。

##土石流発生が流れた範囲
最後に、先ほど保存した画像で白色となった個所のピクセル数を計算します。

import numpy as np
import matplotlib.pyplot as plt
import rasterio

with rasterio.open('images/x-atami_bandnum_02 03 04 08_clusternum_15_20210814-1903_20210805_change_extract_impact.png') as src2:
    data2 = src2.read()

# データのサイズの確認
print('Loaded data'+ str(data2.shape))
plt.imshow(data2[0],cmap='gist_earth')
plt.show()

# 白色ピクセル数
print("White pixel: " + str(np.count_nonzero(data2[0] == 255)))

データが先ほど置き換えた255となっているところのピクセル数を数えます。

Loaded data(3, 160, 200)
White pixel: 245

白色のピクセル数は242個でした。
245 x 10 x 10 = 24,500 m2

つまり、土石流が流れた範囲は、約2万4千平方メートルとなります。

#まとめ

  • 解析結果によると、画像で示された土石流の範囲は約2万4千平方メートル、東京ドーム約半分の広さに影響が及んでいたことが分かりました。
  • 解析画像でも見て取れるように、今回範囲を計算したところは、ほとんどが去年の段階では森林地区として分類されている部分です。そこで起こった2万平方メートル以上の範囲を流れたの土石流が、そのふもとの町に流れ込んだことを考えるといかに被害が大きかったか想像できます。
  • 今回はSentinel2の分解能10mのデータを使用しました。以前使ったランドサットデータは30mだったので、それよりはかなり細かな画像を作ることができ解析もより正確になっていると思います。ただ、それでも1ピクセル10mなのである程度の誤差は出ていると思います。
  • この記事は以前から行っていた機械学習を用いて衛星画像を解析する、IT技術の紹介の目的で書いています。今回は応用してより実社会で起こっている問題を取り上げてみました。今回の記事が、環境問題や防災に役に立つかどうかはわかりませんが何か参考になれば幸いです。

#参考文献および注釈

#最後に
犠牲者の御冥福をお祈り申し上げますとともに、被災された方々に心よりお見舞いを申し上げます。

12
15
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
12
15

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?