記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html

2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

東日本大震災が発生し、今年で9年が経ちます。未曾有の大災害を受けて行政は新規の取り組みに加えて、既存の防災政策の改善を行っています。
その一つに国土交通省によりハザードマップポータルサイトがあります。ポータルサイトは、ハザードマップを利用し、市民の防災意識を高める目的で2007年に作成されたものです。以降、アップデートを重ね、現在では複数の自治体にまたがる複数の災害情報を閲覧することが可能になっています。

ハザードマップを利用すれば、自分の住んでいる場所にどのような災害リスクがあるのかを簡単に確認することができます。
今回は洪水浸水想定区域のハザードマップを開発環境で取得し、Tellus内の衛星データと重ねてみます。さらに、国土数値情報から避難施設データを取得し、取得した画像上に避難場所をプロットしてみましょう。

Tellusの開発環境を使って、オリジナルの防災マップを作成し、自分の住んでいる場所にどのような災害リスクがあるのかを確かめてみましょう。

本記事はTellusの開発環境を用いることを前提としています。
興味のある方はこちらを参考にお申し込みください。

衛星画像と洪水浸水想定区域のタイル画像を取得する

今回の解析で必要となるモジュールをインポートします。

import os, json, requests, math, sys
import numpy as np
from io import BytesIO
import matplotlib.pyplot as plt
from skimage import io
import cv2
import geopandas as gpd
%matplotlib inline

次に、AVNIR-2のタイル画像を取得する関数を準備します。

TOKEN = "自身のトークンを貼り付けてください"

def get_AVNIR_image(zoom, xtile, ytile, opacity=1, r=3, g=2, b=1, rdepth=1, gdepth=1, bdepth=1, preset=None):
    url = "https://gisapi.tellusxdp.com/blend/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format(zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth) 
    if preset is not None:
        url += '&preset='+preset
    headers = {
        "content-type": "application/json",
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

次に、洪水浸水想定区域のタイル画像を取得する関するを準備します。今回は「国管理河川_洪水浸水想定区域（想定最大規模）」のデータを利用しました。

def get_flood_image(z, x, y, params={}):
    """
    浸水域のタイル画像を取得する
    Parameters
    ----------
    z : int
        タイルのズーム率
    x : int
        タイルのX座標
    y : int
        タイルのY座標
    option : dict
        APIパスのオプション(entityId)
    params : dict
        クエリパラメータ
    Returns
    -------
    img: ndarray
        タイル画像
        https://disaportal.gsi.go.jp/hazardmap/copyright/opendata.html#l2shinsuishin_kuni
    """
    url = 'https://disaportaldata.gsi.go.jp/raster/01_flood_l2_shinsuishin_kuni_data/{}/{}/{}.png'.format(z, x, y)
    r = requests.get(url, params=params)
#     print(r.url)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return io.imread(BytesIO(r.content))

上で定義したタイル画像取得関数を利用して、複数のタイル画像をまとめて取得できる関数を作成します。

詳しくは「【ゼロからのTellusの使い方】大きなタイル画像を取得するには？地図タイルを扱う上で知っておきたいTips4選」をご覧ください。

def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, params={}):
    """
    結合されたタイル画像を取得する
    Parameters
    ----------
    get_image_func : function
        タイル画像取得メソッド
        引数はz, x, y, option, params
    z : int
        タイルのズーム率
    topleft_x : int
        左上のタイルのX座標
    topleft_y : int
        左上のタイルのY座標
    size_x : int
         タイルを経度方向につなぎわせる枚数
    size_y : int
        タイルを緯度方向につなぎわせる枚数
    option : dict
        APIパスのオプション(z,x,y除く)
    params : dict
        クエリパラメータ
    Returns
    -------
    combined_image: ndarray
        結合されたタイル画像
    """
    rows = []
    blank = np.zeros((256, 256, 4), dtype=np.uint8)
    for y in range(size_y):
        row = []
        for x in range(size_x):
            try:
                img = get_image_func(z, topleft_x + x, topleft_y + y, params)
            except Exception as e:
                img = blank
            row.append(img)
        rows.append(np.hstack(row))
    return  np.vstack(rows)

今回の解析では東京23区における洪水浸水想定区域のハザードマップを作成します。Tellus OSからタイル座標を取得します。

#以下の数字を変更することで希望のタイル座標に変更できます
z = 12
x = 3636
y = 1611
size_x = 4
size_y = 4

最初に洪水浸水想定区域のタイル画像を取得します。

# 浸水域を取得して描画
combined = get_combined_image(get_flood_image, z, x, y,size_x,size_y)
io.imshow(combined)

続いてAVNIR-2のタイル画像を取得します。取得したタイル画像の範囲は先ほどと同じであることを確認してください。

# X,Yを起点にして、10*10のタイルを取得する
combined_base = get_combined_image(get_AVNIR_image, z, x, y, size_x, size_y)
plt.imshow(combined_base)
plt.show()

AVNIR-2のタイル画像をベースにして、浸水域タイル画像と重ねます。

plt.figure(figsize=(8, 8)) 
blended = cv2.addWeighted(src1=combined_base,alpha=0.8,src2=combined,beta=0.5,gamma=0)
plt.imshow(blended)

上のように別々に処理することもできますが、以下のように処理をまとめる関数を定義すると重なりを確認する際に便利です。

def make_flooded_avnir(zoom, xtile, ytile, xsize, ysize):
    '''AVNIRのタイル画像と地理院の浸水域タイル画像を重ねる'''
    if (zoom < 9) and (zoom > 15): #zoom levelの範囲設定
        sys.exit('zoom level must be lower than 13')
    else:
        z = zoom
        x = xtile
        y = ytile
        flood_img = get_combined_image(get_flood_image, z, x, y, xsize, ysize)
        base_img = get_combined_image(get_AVNIR_image, z, x, y, xsize, ysize)
#      base_imgをベースにして、浸水域を重ねる
        plt.figure(figsize=(8, 8)) 
        blended = cv2.addWeighted(src1=base_img,alpha=0.8,src2=flood_img,beta=0.5,gamma=0)
    return blended

comb_tile_img = make_flooded_avnir(z,x,y,size_x,size_y)
plt.imshow(comb_tile_img)
plt.show()

先ほどと同じ結果が得られていることを確認してください。

このように、Tellusの衛星データと洪水浸水想定区域のデータを簡単に重ね合わせることができました。これだけでも、身の回りのどこが安全そうなのか、ということを把握できるのではないでしょうか？

地図ではなく衛星画像を重ねるメリットとして、最新の周辺情報(建物や道路など)を反映したもので認識できる、ということがあります。
※今回は例のため、過去の衛星画像を重ねています

では、避難施設データを取得し、ハザードマップとどのように重なるのかも確認しましょう。洪水浸水想定区域と避難場所の重なり具合を確かめれば、より安全な避難場所を探すことができます。

タイル画像から緯度経度を算出する関数を定義しましょう。取得したタイル画像から緯度経度を算出できれば、国土数値情報の避難場所データを上手く画像データの上に重ねることができるはずです。

def get_bbox(xtile, ytile, zoom):
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    
    north_east = num2deg(xtile + 1, ytile, zoom)
    south_west = num2deg(xtile, ytile + 1, zoom)

    return (south_west[0], south_west[1], north_east[0], north_east[1])

#保存フォルダを作成
os.makedirs('shelter', exist_ok=True)

国土数値情報から避難施設データを取得・可視化する

東京都の避難場所データは「国土数値情報　避難施設データ」東京都を選択することでダウンロードできます。

同じファイルをダウンロードした場合には、P20-12_13_GMLというシェイプファイルとして保存されます。ダウンロード後、zipを解凍してください。

なお、ダウンロードする際には、アンケートに答える必要があります。必須項目に回答をした後、ファイルをダウンロードしてください。

今回はwork/ディレクトリ内にshelterというフォルダを作成しました。解凍したファイルを全て、このフォルダ内に移動させます。

cpath = os.getcwd()
filepath = cpath+"/shelter/P20-12_13.shp"

では、保存したシェイプファイル(shp)を読み込みましょう。

df_shp = gpd.read_file(filepath,encoding='SHIFT-JIS')
df_shp.head()

タイル座標を緯度経度に変換します。左側のタイル座標が最小経度、下側のタイル画像が最小緯度に相当します。

bbox_nw = get_bbox(x,y,z) #結合された画像の左上のタイル座標の緯度経度
bbox_se = get_bbox(x+size_x-1,y+size_y-1,z) #結合された画像の右下のタイル座標の緯度経度
min_lon = bbox_nw[0]
min_lat = bbox_se[1]

print(bbox_nw)
print(bbox_se)

結合されたタイル画像(comb_tile_img)はピクセル画像の座標を起点に配置されています。そのため原点をゼロとしなければなりません。最小緯度と最小経度をゼロにするために、避難場所の緯度経度をタイル画像の最小緯度と最小経度で引いていきます。

longitudes = []
latitudes = []
for i in range(len(df_shp)):
    longitudes.append(df_shp['経度'][i]-min_lon)
    latitudes.append(df_shp['緯度'][i]-min_lat)

差分をとった後の緯度経度をdiff_longitudes、diff_latitudesとして既存のデータに加えます。

df_shp['diff_longitudes'] = longitudes
df_shp['diff_latitudes'] = latitudes

避難場所のポイントデータを、結合したタイル画像(comb_tile_img)の範囲に収まるようにします。

diff_long_range = bbox_se[2] - bbox_nw[0] # 経度の範囲上限
diff_lat_range = bbox_nw[3] - bbox_se[1] # 緯度の範囲上限
x_img_size = comb_tile_img.shape[0] #画像のサイズを求める（x軸方向）
y_img_size = comb_tile_img.shape[1] #画像のサイズを求める（y軸方向）

1ピクセルあたりの緯度経度変化を求めます。

onepix_degree_long = diff_long_range/x_img_size
onepix_degree_lat = diff_lat_range/y_img_size
#1ピクセルあたりの緯度経度変化
print(onepix_degree_lat)
#0.0002788934872981
print(onepix_degree_long)
#0.00034332275390625

取得したタイル画像に含まれる範囲の避難場所のみを抽出します。

df_ext = df_shp[(df_shp.diff_longitudes >= 0) & (df_shp.diff_latitudes >= 0) & \
                (df_shp.diff_longitudes <= diff_long_range) & (df_shp.diff_latitudes <= diff_lat_range)]

df_ext.head()

範囲を絞った避難場所のプロットしてみましょう。

plt.scatter(df_ext.diff_longitudes, df_ext.diff_latitudes)
plt.show()

避難場所の緯度経度の値をピクセル座標の値に変換します。これにより、衛星画像のデータに避難場所を重ねることが可能になります。

pixel_x = df_ext.diff_longitudes/onepix_degree_long
pixel_y = df_ext.diff_latitudes/onepix_degree_lat
pixel_x = np.array(pixel_x,dtype = 'int16')
pixel_y = np.array(pixel_y,dtype = 'int16')

算出したピクセル座標をデータフレームに加えます。

df_ext.insert(df_ext.shape[1], "pixel_x", pixel_x, True) 
df_ext.insert(df_ext.shape[1], "pixel_y", pixel_y, True) 
df_ext.head()

再度描画してみましょう。緯度経度がピクセル座標に変換されていることが確認できます。

plt.figure(figsize=(7, 7))
plt.xlim(0,x_img_size)
plt.ylim(0,y_img_size)
plt.scatter(df_ext.pixel_x, df_ext.pixel_y)
plt.show

最後に取得済みのタイル画像と避難場所の位置データを重ねましょう。

fig, ax = plt.subplots(figsize=(8, 8))
blended = np.flipud(comb_tile_img)#画像の上下を反転する（プロット範囲の指定により、配列が反転してしまうため）
ax.imshow(blended)
plt.xlim(0,x_img_size)
plt.ylim(0,y_img_size)
plt.scatter(df_ext.pixel_x, df_ext.pixel_y, s=5, alpha=0.8)
plt.show

一連の流れを関数で定義する

今回は東京都のシェイプファイルを利用していますので、他の都道府県における避難場所示したい場合は関数のshp_fileを変える必要があります。

def flood_shelter(shp_file, z, x, y, size_x=1, size_y=1,psize=5,palpha=0.8):
    '''衛星画像と洪水浸水想定区域、さらに避難場所を重ねた画像を出力する関数
    -------------------------------
    Parameters
    -------------------------------
    shp_file: 参照するシェイプファイル
    z: zoom level(9-14)
    x: タイルのX座標
    y: タイルのY座標
    size_x: 取得するタイルの数(x軸方向)。デフォルト1
    size_y: 取得するタイルの数(y軸方向)。デフォルト1
    psize: ポイントのサイズ。デフォルト5
    palpha: ポイントの不透明度デフォルト0.8
    '''
    comb_tile = make_flooded_avnir(z,x,y,size_x,size_y) #衛星と洪水浸水想定区域の結合
    df_shp = gpd.read_file(filepath,encoding='SHIFT-JIS') #シェイプファイルの読み込み
    bbox_nw = get_bbox(x,y,z) #結合された画像の左上のタイル座標の緯度経度
    bbox_se = get_bbox(x+size_x-1,y+size_y-1,z) #結合された画像の右下のタイル座標の緯度経度
    min_lon = bbox_nw[0] #最小経度
    min_lat = bbox_se[1] #最小緯度
    longitudes = []
    latitudes = []
    for i in range(len(df_shp)):
        longitudes.append(df_shp['経度'][i]-min_lon)
        latitudes.append(df_shp['緯度'][i]-min_lat)
    
    
    df_shp['diff_longitudes'] = longitudes
    df_shp['diff_latitudes'] = latitudes
    diff_long_range = bbox_se[2] - bbox_nw[0] # 経度の範囲上限
    diff_lat_range = bbox_nw[3] - bbox_se[1] # 緯度の範囲上限
    x_img_size = comb_tile_img.shape[0] #画像のサイズを求める（x軸方向）
    y_img_size = comb_tile_img.shape[1] #画像のサイズを求める（y軸方向）
    
    #1ピクセルあたりの緯度経度変化を求める
    onepix_degree_long = diff_long_range/x_img_size
    onepix_degree_lat = diff_lat_range/y_img_size
    
    #避難場所を取得する範囲を制限する
    df_ext = df_shp[(df_shp.diff_longitudes >= 0) & (df_shp.diff_latitudes >= 0) & \
                (df_shp.diff_longitudes <= diff_long_range) & (df_shp.diff_latitudes <= diff_lat_range)]
    
    #1ピクセルあたりの座標値を求める
    pixel_x = df_ext.diff_longitudes/onepix_degree_long
    pixel_y = df_ext.diff_latitudes/onepix_degree_lat
    pixel_x = np.array(pixel_x,dtype = 'int16')
    pixel_y = np.array(pixel_y,dtype = 'int16')
    
    #新しい変数をデータフレームの最後尾に加える
    df_ext.insert(df_ext.shape[1], "pixel_x", pixel_x, True) 
    df_ext.insert(df_ext.shape[1], "pixel_y", pixel_y, True) 

    #画像の描画
    fig, ax = plt.subplots(figsize=(8, 8))
    blended = np.flipud(comb_tile)#画像の上下を反転する（プロット範囲の指定により、配列が反転してしまうため）
    ax.imshow(blended)
    plt.xlim(0,x_img_size)
    plt.ylim(0,y_img_size)
    plt.scatter(df_ext.pixel_x, df_ext.pixel_y, s=psize, alpha=palpha)
    plt.show()

実行すると、先ほどと同じ画像を取得できることを確認してください。

flood_shelter(filepath,z,x,y,size_x,size_y)

ズームレベルとタイル座標を変化させると、より詳細な場所を確認することができます。

z,x,y = (13,7276,3224)
flood_shelter(filepath,z,x,y,size_x,size_y,psize = 15)

以上で、開発環境を利用したオリジナルハザードマップの作成は完了しました。

ハザードマップポータルでは、洪水浸水想定区域に他にも土砂災害や津波などについてのデータも提供しています。

本記事と同じ手順でAPIを利用し、他のデータも取得することが可能です。開発環境でマイハザードマップ作成に取り組んでみてください。

番外編：ローカル環境でマイマップを作成する

ここからは番外編になります。

番外編では、衛星画像と洪水浸水想定区域のpng画像をgeotiffに変換し、QGISを始めとしたGISアプリケーションでも利用できるようにしてみましょう。
png画像からgeotiffを作成するための関数を定義します。

from math import log,tan,radians,cos,pi,floor,degrees,atan,sinh
import glob, subprocess
import gdal # subprocess, gdalをインストールしていなければ、pip install gdalとpip install subprocessをterminal上で実行してください

# 以下のコードはhttps://jimmyutterstrom.com/blog/2019/06/05/map-tiles-to-geotiff/を参考にしています。

最初にタイル座標を経度データに変換するための関数を定義します。

def x_to_lon_edges(x, z):
    tile_count = pow(2, z)
    unit = 360 / tile_count
    lon1 = -180 + x * unit
    lon2 = lon1 + unit
    return(lon1, lon2)

緯度データにも同様の処理をします。この際にメルカトル図法に合わせて緯度を変換していきます。

def mercatorToLat(mercatorY):
    return(degrees(atan(sinh(mercatorY))))


def y_to_lat_edges(y, z):
    tile_count = pow(2, z)
    unit = 1 / tile_count
    relative_y1 = y * unit
    relative_y2 = relative_y1 + unit
    lat1 = mercatorToLat(pi * (1 - 2 * relative_y1))
    lat2 = mercatorToLat(pi * (1 - 2 * relative_y2))
    return(lat1, lat2)

上で定義した関数を用いて、取得した座標を配列として返します。

これはgdal.Translateでtiff画像に空間参照系であるepsg4326を当てはめるために必要なプロセスになります。空間参照系を画像のデータに割り当てることにより、画像データは正確に地理座標上に重ねられることになります。空間参照系について詳しくはこちらをご覧ください。

def tile_edges(x, y, z):
    lat1, lat2 = y_to_lat_edges(y, z)
    lon1, lon2 = x_to_lon_edges(x, z)
    return[lon1, lat1, lon2, lat2]

最後にタイル画像をtiffに変換する関数を定義します。

def georeference_raster_tile(x, y, z, path):
    bounds = tile_edges(x, y, z)
    filename, extension = os.path.splitext(path)
    gdal.Translate(filename + '.tif',
                   path,
                   outputSRS='EPSG:4326',
                   outputBounds=bounds)

それでは関数を実行してみましょう。最初に抽出するタイル画像の範囲を指定します。

zoom, x_min, x_max, y_min, y_max = (12, 3636, 3639, 1611, 1614)

def make_single_flooded_avnir(zoom, xtile, ytile):
    '''AVNIRのタイル画像と地理院の浸水域タイル画像を重ねる'''
    if zoom > 15:
        sys.exit('zoom level must be lower than 14')
    else:
        z = zoom
        x = xtile
        y = ytile
        flood_img = get_combined_image(get_flood_image,z, x, y)
        base_img = get_AVNIR_image(z, x, y)
#      base_imgをベースにして、浸水域を重ねる
        plt.figure(figsize=(8, 8)) 
        blended = cv2.addWeighted(src1=base_img,alpha=0.8,src2=flood_img,beta=0.5,gamma=0)
        return blended

#保存フォルダを作成
os.makedirs('avnir_imgs', exist_ok=True)

for x in range(x_min, x_max + 1):
    for y in range(y_min, y_max + 1):
        print(f"{x},{y}")
        img = make_single_flooded_avnir(zoom,x,y)
        plt.imsave(os.getcwd()+'/avnir_imgs/'+str(x)+'_'+str(y)+'.png',img)
        png_path = os.getcwd()+'/avnir_imgs/'+str(x)+'_'+str(y)+'.png'
        georeference_raster_tile(x, y, zoom, png_path)

取得し、変換したgeotiff画像を一枚の画像に結合します。

def merge_tiles(input_pattern, output_path):
    merge_command = ['gdal_merge.py', '-o', output_path]

    for name in glob.glob(input_pattern):
        merge_command.append(name)

    subprocess.call(merge_command)

tiff_dir = os.getcwd()+'/avnir_imgs/'

merge_tiles(tiff_dir + '/*.tif', tiff_dir + '/merged.tif')

結合した画像を確認しましょう。

import rasterio
from rasterio.plot import show

src = rasterio.open(tiff_dir+'merged.tif')
plt.figure(figsize = (8,8))
show(src.read(), transform=src.transform)

冒頭で取得した、衛星画像と洪水浸水想定区域を重ねた画像と同じ画像を取得できるかと思います。見た目は同じですが、画像をダウンロードしてみると、画像の形式が異なることが分かるかと思います。

実際にQGISで画像を確認してみましょう。確認方法については例えばこちらの記事をご覧ください。

以下では操作手順をスクショで紹介します。

以上が番外編、png画像からgeotiff画像を作成してみようでした。

まとめ

今回の記事では3/11という日に、安全について考えるきっかけとして、浸水予想域の情報や避難所の情報を取得する方法をご紹介しました。

流れをおさらいすると、
・Tellusに搭載された衛星画像を取得する方法
・Tellus環境外の画像を取得する方法
・データ同士を重ね合わせる方法
・シェイプファイルを重ね合わせる方法
を紹介しました。

今回の記事を通して、身の回りの避難所の情報含め、公開されているハザードマップなどを確認し、発災時に適切に行動できるように備えて頂ければ幸いです。