推計震度分布図、ほしいですか？

タイトルめが突然失礼しました。
始めまして。地震情報アプリ「Zero Quake」を開発しております、べにだてです。
皆さんが欲しくてたまらないであろう推計震度分布図。震度の分布が視覚的に表現され、被害がわかりやすい推計震度分布図。今回はそんな推計震度分布のお話です。

推計震度分布図ってなんぞや

さっきからくどいほどに登場する、絶妙に親しみにくい「推計震度分布図」という言葉。

推計震度分布図＝観測点で観測された震度を、面的に補完した図

つまり、観測点の無いところも含めて、その地震による揺れを推定したものです。日本全土を250m四方¹の四角形（メッシュ）に分割し、四角形毎に震度が推計されています。気象庁からは、震度４以上の地域を推計した情報が提供されます。

↓実際の推計震度分布図（出典：気象庁ホームページ）

震度の推計には、震度観測点で計測した震度に加え、その土地の揺れやすさが反映されます。
この推計震度分布図は、主に最大震度５弱以上の地震の際に発表されます。しかし、近年、すべての有感地震に推計震度分布図を提供しようという案も出ているようです。（気象庁資料で見た気がします）

推計震度分布図についてさらに詳しく知りたい方は、気象庁ホームぺージを確認してください。

データを探して三千里

地震関係の情報を取得する、というと、気象庁防災情報XMLが真っ先に浮かびます。しかし、ここでは推計震度分布図は提供されません。推計震度分布図を提供する無料のAPIもありません。
八方ふさがり...いえ、違います。気象庁ホームページから取得すればいいのです！

気象庁ホームページの情報は、出典を記載すれば、自由に複製や公衆送信ができます。
詳しくは利用規約を参照。

次の章では、気象庁ホームぺージから推計震度分布図を取得する方法を紹介します。

気象庁HPからデータを取ろう

0⃣ 前提知識　データ形式

気象庁ホームページでは、推計震度分布が画像として提供されています。
画像と言っても、1枚の画像ではなく、1次メッシュごとに分割され、800×800pxの透過png画像となっています。
データは「Webメルカトル図法」で投影されており、地理院地図やOpenStreetMap等と重ねることができます。
6つの震度階級が色で示されており、震度4以下の地域は、透明となっています。

1⃣ リストを取得

最初のステップは、「リストの取得」です。
過去に発表された推計震度分布のリストが、以下のURLで提供されています。驚くべきことに、過去の情報が すべて （現時点で270件程）記載されています。

https://www.jma.go.jp/bosai/estimated_intensity_map/data/list.json

このリストからは、震源要素、最大震度、情報が提供されている範囲の１次メッシュコード等の情報が得られます。

list.json

[{
		"hypo": { //震源要素など
			"it": "2023-06-11T19:00:06",	//発表時刻？
			"at": "2023-06-11T18:55:00",	//地震検知日時
			"lat": 42.54,					//震源の緯度
			"lon": 141.95,					//震源の経度
			"dep": 140,						//震源の深さ
			"mag": 6.2,						//マグニチュード
			"epi": "浦河沖",				//震央地名
			"kun": 0,
			"maxi": 4.5						//最大震度
		},
		"comment": "震度４の地域では、座りの悪い置物が倒れるなどしている可能性があります。",
		"url": "202306111855_192",			//【Ⓐ】 URLの生成に使用
		"rank_cnt": {"i9": 0, "i8": 0, "i7": 0, "i6": 0, "i5": 163, "i4": 56573, "i3": 497901, "i2": 682224, "i1": 535730, "i0": 26929},
		"bounds": [				//データの範囲
			[40.45, 140.94],	//南西の端
			[43.4, 145.13]		//北東の端
		],
        //↓ 【Ⓑ】 データを含むメッシュコード（URLの生成に使用）
		"mesh_num": ["6041", "6141", "6240", "6241", "6341", "6342", "6343","6441", "6442", "6443", "6444","6545"
		],
		"datum": 2
},{
    //省略...
}]

【Ⓐ】、【Ⓑ】で示した情報は、以降のステップで使用します。

参考情報

json[n].hypo.itは、2000-01-01T00:00:00（おそらくデータなしの意）となっているケースもあります。
json[n].hypo.kunの値は、現時点で常に 0 です。

2⃣ 画像のURLを生成

画像はどこだ！リストにURLなんてなかったじゃないか！と言いたくなった方、しばらくお待ちください。今から、肝心の画像のURLを生成します。

まずはこちらをご覧ください。

https://www.jma.go.jp/bosai/estimated_intensity_map/data/【Ⓐ】/【Ⓑ】.png

ここに、リストから取得した情報を以下のように当てはめます。※上の章の記号と対応します。

	当てはめる値	説明	例
【Ⓐ】	`json[n].url`	情報のID	"202306111855_192"
【Ⓑ】	`json[n].mesh_num[m]`	1次メッシュコード	"6041"

https://www.jma.go.jp/bosai/estimated_intensity_map/data/202305051442_390/5636.png

おめでとうございます、やっとURLが完成しました。「mesh_num」に含まれるメッシュコードの数だけ（今回の例では15回）、URLを生成しておきましょう。

←取得できる画像の例（出典：気象庁ホームページ）

3⃣ 表示

画像が取得できたら、地図上に表示します。

準備するもの

LeafletやMapLibre GL JSなど、お好きな地図ライブラリ
画像のURLとメッシュコード × n個

⑴ 緯度・経度の取得

画像を地図上に描画するには、画像のほかに、画像の四方の座標が必要です。今回は、メッシュコードからこの座標を計算します。

前提知識１次メッシュコードの規則

前半２桁…南端の緯度 × 1.5
後半２桁…西端の経度 - 100

東端の経度 = 後半２桁 + 100 + 1;
西端の経度 = 後半２桁 + 100;
南端の緯度 = 前半２桁 / 1.5;
北端の緯度 = 前半２桁 / 1.5 + 2/3;

上のように、メッシュコードには規則があるため、コードから緯度経度を求めることができます。

mesh.js

var MeshCode = "6041";
var lat = Number(MeshCode.substring(0, 2)) / 1.5;
var lng = Number(MeshCode.substring(2, 4)) + 100;
var lat2 = lat + 2 / 3;
var lng2 = lng + 1;

おっと、私が未だにletではなくvarを使っていることがバレてしまいました。

⑵ 緯度・経度の取得

四方の座標がわかったら、あとはライブラリに任せて描画するだけです。
ここでは、LeafletとMapLibreでの描画の例を挙げます。

MapLibreでの描画例

MapLibre_EstimatedIntensityMap.js

var image_url = "https://www.jma.go.jp/bosai/estimated_intensity_map/data/202203162336_289/5640.png";
var map = new maplibregl.Map({ container: "map" });

map.addSource("estimated_intensity_map", {
    type: "image",
    url: image_url,
    coordinates: [  //先ほど計算した四方の座標
        [lng, lat2],
        [lng2, lat2],
        [lng2, lat],
        [lng, lat],
    ],
});

map.addLayer({
    id: "estimated_intensity_map_layer",
    type: "raster",
    source: "estimated_intensity_map",
    paint: {
        "raster-resampling": "nearest", //拡大時にぼかさない
    },
});

MapLibreでは、「raster-resampling」オプションを「"nearest"」に設定することをお勧めします。こうすることで、拡大時に、メッシュの境界がぼかされることがなくなります。²これは、CSSのimage-renderingプロパティを "pixelated"に設定するのと似ています。

Leafletでの描画例

Leaflet_EstimatedIntensityMap.js

var image_url = "https://www.jma.go.jp/bosai/estimated_intensity_map/data/202203162336_289/5640.png";
var map = L.map("map");

L.imageOverlay(
    image_url,
    [[lat2, lng], [lat, lng2]] //先ほど計算した四方の座標
).addTo(map);

データの利用上の注意点

１ピクセル = １メッシュではない

画像1枚 = 1次メッシュの範囲を 800×800px で表現

	メッシュ解像度	画像1枚のメッシュ数	1メッシュあたり
以前	1km	80×80	10px
現在	250m	320×320	2.5px ←⁉

上の表を見てわかる通り、現在提供されている250mメッシュの推計震度分布図は、800×800pxで表現しようとすると、1メッシュあたりのピクセル数が小数になってしまいます。すると、想像がつくかもしれませんが...

上の画像のように、メッシュの境目がぼやけてしまうわけです。こればかりは、気象庁側のシステムの問題なので、どうしようもありません。
相当ズームしない限り弊害は生まれないはずですが...

しかし、先日、これに対する対策法を編み出しました。続報をお待ちください。

1kmメッシュ…1次メッシュの80分の1
(1次メッシュを8分割した2次メッシュを10分割した3次メッシュと同じ）

250mメッシュ…1次メッシュの320分の1
(1次メッシュを8分割した2次メッシュを10分割した3次メッシュを4分割したもの)

（ややこしい）

2023年1月以前は1km四方のメッシュで提供されていました。 ↩
現在は、元のデータの時点で境界がぼやけています。詳しくはデータの利用の章を参照。 ↩