令和6年能登半島地震につきまして、被災された皆さまに心からお見舞い申し上げます。

災害は平時の備えにより被害を大きく減らすことができます。
今回は自分の住むつくば市について「地盤の強さ」すなわち「揺れやすさ」をマップにして地域の地震の影響に対する理解を高めたいと思います。

今回は防災科学技術研究所が公開している表層地盤増幅率を使います。同データは日本の全国の地域についてまとめられていますので、つくば市だけでなく任意の市町村についてマップをつくることができます。なお、お住まいの市町村についても、自治体が大学や国の研究機関と連携してハザードマップとして整備している場合が多いかと思います。

さて、表層地盤増幅率という言葉について、耳慣れない方が多いかと思いますのでWikiで調べてみましょう。

簡単に言えば土地の揺れやすさを数字にしたものです。
地域ごとの表層地盤増幅率は防災科研の公開するデータから入手できます。

つくば市全域については、5339、5340、5439、5440の4つのファイルに含まれるデータが該当しそうですので、全てDLします。

データはかなり広い範囲の地域を含みますので、つくば市にあわせて経度緯度を絞ります。大体、経度は139.9~140.3の間、緯度は35.9~36.3の間くらいでしょうか。

まずは表層地盤データとつくば市の概形を重ねるところまでをPythonで実装してみます。いつも通りgeopandasをインストールし、必要なライブラリをインポートします。

!pip install --upgrade geopandas --q

import pandas as pd
import numpy as  np
import matplotlib.pyplot as plt
import geopandas as gpd
import folium
from branca.element import Figure

使用するファイルを読み込みます。

dir_path = 'データの置き場所'

gdf_area = gpd.read_file(dir_path+'/N03-23_08_230101.shp',encoding='shift-jis')
gdf_area = gdf_area[gdf_area['N03_004']=='つくば市']

gfiles = [dir_path+'/表層地盤/5339/Z-V4-JAPAN-AMP-VS400_M250-SHAPE-5339.shp',
          dir_path+'/表層地盤/5340/Z-V4-JAPAN-AMP-VS400_M250-SHAPE-5340.shp',
          dir_path+'/表層地盤/5439/Z-V4-JAPAN-AMP-VS400_M250-SHAPE-5439.shp',
          dir_path+'/表層地盤/5440/Z-V4-JAPAN-AMP-VS400_M250-SHAPE-5440.shp']
gdf =  gpd.GeoDataFrame(pd.concat([gpd.read_file(files, encoding='Shift-JIS') for files in gfiles]))

データの範囲を絞り込んで、重畳して表示します。

maxx = 140.22
minx = 139.95
maxy = 36.28
miny = 35.90

gdf = gdf[(gdf['geometry'].bounds.minx >= minx)&
　　　　　(gdf['geometry'].bounds.maxx <= maxx)&
　　　　　(gdf['geometry'].bounds.miny >= miny)&
　　　　　(gdf['geometry'].bounds.maxy <= maxy)]

gdf["geometry"] = gpd.GeoSeries(gdf["geometry"]).simplify(tolerance=0.001)
gdf= gdf.reset_index(drop=True)

fig = Figure(width=300, height=300)
map = folium.Map(location=[36, 140] ,zoom_start=9)

area_style_function = lambda x: {'color' :  'blue','opacity' : 0.70,'weight' : 1,}
folium.GeoJson(gdf_area,style_function = area_style_function).add_to(map)

box_style_function = lambda x: {'color' :  'red','opacity' : 0.50,'weight' : 3,}
folium.GeoJson(gdf['geometry'].unary_union, style_function = box_style_function).add_to(map)
map

表層地盤データと市の概形のデータが大体重なりました。
最後に、概形に合わせて表層地盤データを切り取り、表層地盤増幅率に合わせて色分けをします。

まずは市の概形データから境界ポリゴンデータを作成します。

boundary_polygon = gdf_area.unary_union

数字に合わせてカラーマップを作ります。

import matplotlib.cm as cm

v= np.arange(0.5, 2.5, 0.25)
c = cm.viridis(range(0,250,int(250/len(v))))
fig, ax = plt.subplots(figsize = (10,6),dpi=100)
ax.imshow([c])

いつも通り、’色’というカラムを作り、それぞれのメッシュに対応するHTMLカラーの値を格納します。

for i in range(0,len(gdf)):
  r = float(gdf.iloc[i]['ARV'])

  if r == 0:
    gdf.loc[i, '色'] = '#000000'
    gdf.loc[i, '色'] = '#0000'

  else:
    id = sum(v >= r)
    if id == 0:
      r=c[-1]
    else:
      r = c[-id]
    r = tuple((r[0:3]*255).astype(int))
    r = rgb_to_hex(r)
    gdf.loc[i, '色'] = '#'+str(r)

マップを作成します。

fig = Figure(width=900, height=1200)
map = folium.Map(location=[36.1, 140.1], tiles='cartodbdark_matter',zoom_start=12,zoom_snap=0.2)

fig.add_child(map)

area_style_function = lambda x: {'color' :  'blue','opacity' : 0.00,'weight' : 0,}
folium.GeoJson(gdf_area,style_function = area_style_function).add_to(map)

sf = lambda x:{"fillColor": x['properties']['色'],
               "color": '#000000',
               "weight": 0.5,
               "fillOpacity": 0.9}
folium.GeoJson(gdf[gdf["geometry"].intersects(boundary_polygon)],style_function=sf).add_to(map)

市の中心部の表層地盤増幅率はおおむね2.0程度であり、防災技術研究所の指標でみれば「揺れやすい」地域に分類されます。市の南東側の小貝川付近は特に揺れやすい地区であり、注意が必要でしょう。市の中心部にも2.25~2.5の地域があります。反対に、筑波山近くになるとかなり丈夫な地盤となっているようです。

最後にこれまでにも使った河川や鉄道、道路のデータについても重畳してみてみます。

作成したマップをHTMLで保存しておきます。

map.save(dir_path+'map.html')