防災アプリ開発 Advent Calendar 2023 12日目
経緯
地震情報アプリ「Zero Quake」を制作しているべにだてと申します。
現状、無料で推計震度分布図を取得する方法は、ほぼ気象庁ホームページ一択です。しかし、ここで公開されている情報は、扱いづらい画像形式。今回は、この気象庁ホームページから得られる画像を、扱いやすいGeoJSON形式に変換するべく、四苦八苦したお話です。
ちなみに、無料でベクター形式の推計震度分布図を取得する方法のご紹介はおそらく世界初です（しょうもない）
気象庁HPのコンテンツは、利用規約により、出典を記載すれば、ある程度自由に利用できます。この記事の内容を実践する際は、ご自身で規約をよく確認してください。
当記事の手法は、当方で十分確認してはおりますが、誤ったデータを出力する可能性があります。自己責任でご利用ください。
気象庁HPから推計震度分布図を取得する
気象庁HPより ／ 推計震度分布図の一例
今回は、気象庁ホームページをスクレイピングして、推計震度分布図を取得します。
推計震度分布図は、１次メッシュ1毎に分割されており、取得方法は複雑です。詳細な取得方法は、「第三編 推計震度分布への道」を参照してください。
今回は、画像の取得の手順を省いて、変換の手順を紹介します。
GeoJsonにしてしまえ！！
画像データをGeoJSON形式に変換すれば、何かと扱いやすくなります。しかし、この変換、一筋縄にはいきません。
変換の壁
気象庁ホームページでは、推計震度分布図は画像として提供されます。言わずもがな、画像というのは、プログラムからは扱いづらいものです。詳細を説明します。
【悲劇】１ピクセル = １メッシュではない
気象庁HPの推計震度分布図は、1次メッシュ1ごとに分割され、800×800px のラスター画像で表現されています。そして、画像１枚に注目すると...
| メッシュ間隔 | 解像度 | メッシュ１個あたり | |
|---|---|---|---|
| 以前 | 1km | 1px=100m | 10px |
| 現在 | 250m | 1px=100m | 2.5px ←⁉ |
上の表を見てわかる通り、現在提供されている250mメッシュ1の推計震度分布図は、800×800pxで表現しようとすると、メッシュ1つあたりのピクセル数が小数になってしまいます。すると、想像がつくかもしれませんが...
上の画像のように、アンチエイリアスの作用で、メッシュの境目がぼやけてしまうわけです。明らかにプログラムから処理しにくそうですね。こればかりは、気象庁側のシステムの問題なので、どうしようもありません。この特徴が、様々な処理に影響してしまいます....
どうやって変換する？
まず、気象庁HPの画像と元データを重ねてみましょう。重ねた模式図が以下です。
１ピクセル=１メッシュではない ため、アンチエイリアスにより、図中の２番、７番のように、震度階級に属さない、「微妙な色」が発生してしまいます。
しかし、心配は不要。メッシュ１つ分（250m）は2.5ピクセルで表されます。2.5>2であるため、必ず、上の図で言う、0番,1番,3番,6番,8番...のような、元データの色を完全に反映しているピクセルがあるわけです。
つまり、このピクセルのみをかいつまんで取得していけば...
上のように、元データが復元できるわけです！！
そして、この「元の色を反映しているピクセル」の分布には、法則があります。これまでは、x軸のみを見て説明していましたが、実際には、y軸にも同じことが言えます。
そして、その法則が、以下の通り。
（１時間ほどかけて確認しました）
for文を使い、工夫してピクセルの色を取得して、GeoJSONを生成すれば...なんだかうまくいきそうですね！
変換の実装
実際に実装していきます。使用言語は、ブラウザで動作するJavascript（Vanilla JS）です。今回は、色を解析して、震度を割り出す機能も実装してみます。
Step.1 まずは取得
var url = "202305051442_390"; //情報のID
var MeshCode = "5637"; //画像の１次メッシュコード
var EstimateShindoData = [];
//画像の四方の範囲
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; //画像東端の経度
var ESMap_canvas = document.createElement("canvas");
ESMap_canvas.width = ESMap_canvas.height = 800;
var context = ESMap_canvas.getContext("2d");
var img = new Image();
img.onload = function () {
context.clearRect(0, 0, 800, 800);
context.drawImage(img, 0, 0, 800, 800); //画像を描画
/*ここに次の章で変換処理を書いていくよ！*/
};
//推計震度分布図のURLを生成
img.src = "https://www.jma.go.jp/bosai/estimated_intensity_map/data/" + url + "/" + MeshCode + ".png";
まずは、気象庁HPから、推計震度分布図の画像を取得します。取得した図は、処理に備えて、canvasに描画しておきます。
今回は、サンプルプログラムですので、画像の取得は簡易的に実装しています。取得についての詳細は、以前投稿した「第三編 推計震度分布への道」を参照してください。
また、画像の東西南北の端の座標も求めておきます。
Step.2 GeoJSONに加工
//...略
img.onload = function () {
context.clearRect(0, 0, 800, 800);
context.drawImage(img, 0, 0, 800, 800);
var y = 1;
var imgData = context.getImageData(0, 0, 800, 800).data;
for (let i = 0; i < 320; i++) {
var x = 1;
for (let j = 0; j < 320; j++) {
var r = imgData[(y * 800 + x) * 4];
var g = imgData[(y * 800 + x) * 4 + 1];
var b = imgData[(y * 800 + x) * 4 + 2];
var a = imgData[(y * 800 + x) * 4 + 3];
var int;
if (a > 50) {
if (Math.abs(r - 250) < 16 && Math.abs(g - 230) < 16 && Math.abs(b - 150) < 16) int = "4";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 230) < 16 && Math.abs(b - 0) < 16) int = "5-";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 153) < 16 && Math.abs(b - 0) < 16) int = "5+";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 40) < 16 && Math.abs(b - 0) < 16) int = "6-";
else if (Math.abs(r - 165) < 16 && Math.abs(g - 0) < 16 && Math.abs(b - 33) < 16) int = "6+";
else if (Math.abs(r - 180) < 16 && Math.abs(g - 0) < 16 && Math.abs(b - 104) < 16) int = "7";
var North = lat2 + ((lat - lat2) / 320) * i;
var South = lat2 + ((lat - lat2) / 320) * (i + 1);
var West = lng + ((lng2 - lng) / 320) * j;
var East = lng + ((lng2 - lng) / 320) * (j + 1);
EstimateShindoData.push({
type: "Feature",
properties: { rgb: [r, g, b], Intensity: int },
geometry: {
type: "Polygon",
coordinates: [
[
[West, North],
[East, North],
[East, South],
[West, South],
[West, North],
],
],
},
});
}
if (j % 2 == 0) x += 3;
else x += 2;
}
if (i % 2 == 0) y += 2;
else y += 3;
}
var geojson = {
type: "FeatureCollection",
features: EstimateShindoData,
};
};
//略...
各ピクセルの色を取得し、GeoJSONを生成するコードです。各部を順に説明します。
Step.2-1 色を取得
var y = 1;
var imgData = context.getImageData(0, 0, 800, 800).data;
for (let i = 0; i < 320; i++) {
var x = 1;
for (let j = 0; j < 320; j++) {
//変数x=色を取得するX座標 変数y=色を取得するY座標
/*...略...*/
if (j % 2 == 0) x += 3;
else x += 2;
}
if (i % 2 == 0) y += 2;
else y += 3;
}
まず、二重のfor文で、左から右へ、上から下へ、順にピクセルを参照します。先ほどの説明の通り、かいつまんで参照します。そのため、規則に沿って、変数を交互に２増加させたり、３増加させたりしています。
Step.2-2 色を震度に変換
var r = imgData[(y * 800 + x) * 4]; //ピクセルの色 レッド
var g = imgData[(y * 800 + x) * 4 + 1]; //ピクセルの色 グリーン
var b = imgData[(y * 800 + x) * 4 + 2]; //ピクセルの色 ブルー
var a = imgData[(y * 800 + x) * 4 + 3]; //ピクセルの色 透明度
RGBA形式で、参照したピクセルの色を取得します。
var int;
if (a > 50) { //ピクセルの色が透明ではなかったら
if (Math.abs(r - 250) < 16 && Math.abs(g - 230) < 16 && Math.abs(b - 150) < 16) int = "4";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 230) < 16 && Math.abs(b - 0) < 16) int = "5-";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 153) < 16 && Math.abs(b - 0) < 16) int = "5+";
else if (Math.abs(r - 255) < 16 && Math.abs(g - 40) < 16 && Math.abs(b - 0) < 16) int = "6-";
else if (Math.abs(r - 165) < 16 && Math.abs(g - 0) < 16 && Math.abs(b - 33) < 16) int = "6+";
else if (Math.abs(r - 180) < 16 && Math.abs(g - 0) < 16 && Math.abs(b - 104) < 16) int = "7";
//略...
}
ピクセルの色と各階級の色と比較して、震度階級を判断します。
画像の圧縮などで、色が少しばかり変わってしまうことがあるため、RGBの各値を比較し、それぞれ差分16未満を等しい値とみなします。
（例）Math.abs(r - 【定数】) < 16 … 「変数r の値が【定数】とほぼ等しいなら」という意味
Step.2-3 GeoJSONポリゴンを作成
ここまでで、メッシュの色と震度を取得できました。あとは、これをGeoJSONに変換します。
var North = lat2 + ((lat - lat2) / 320) * i; //250mメッシュの北端の緯度
var South = lat2 + ((lat - lat2) / 320) * (i + 1); //250mメッシュの南端の緯度
var West = lng + ((lng2 - lng) / 320) * j; //250mメッシュの西端の経度
var East = lng + ((lng2 - lng) / 320) * (j + 1); //250mメッシュの右端の経度
計算によって、該当する250メッシュ1の東西南北の端の座標を求めます。これが、GeoJSON Featureの形のもととなります。
1次メッシュ1を縦横に320分割し、現在のループ回数などから座標を計算しています。
//GeoJSON Feature（ポリゴン）を作成、配列に追加
EstimateShindoData.push({
type: "Feature",
properties: { rgb: [r, g, b], Intensity: int }, //色・震度をpropertiesに保存
geometry: {
type: "Polygon",
coordinates: [ //メッシュのジオメトリ
[
[West, North],
[East, North],
[East, South],
[West, South],
[West, North],
],
],
},
});
先ほど計算した端の座標をもとに、５点の座標を定義し、多角形のポリゴンを作成します。propertiesに、取得したRGBカラーや震度を格納します。
④ ひとつのGeoJSONにまとめる
var geojson = {
type: "FeatureCollection",
features: EstimateShindoData,
}
先ほどのステップで、メッシュ毎にfeatureを作成しました。これをFeatureCollectionにまとめて、GeoJSONの完成です！あとは、お好きなライブラリで描画するなり、ファイルに出力するなり、ご自由にどうぞ。
ただし、今回処理したのは、画像１枚分です。気象庁HPの画像は１次メッシュ1ごとに分割されており、たいていの場合数枚に分かれています。ループして取得・解析する処理などは、ご自身で実装してください。
他の言語で実装したい方へ
私Javascript読めない...そんなあなたへ！他の言語で実装したい方のために、実装の大枠を簡単に示しておきます。
lat = /*メッシュコードの前半２桁*/ / 1.5; //画像西端の経度
lng = /*メッシュコードの後半２桁*/ + 100; //画像北端の緯度
lat2 = lat + 2 / 3; //画像南端の緯度
lng2 = lng + 1; //画像東端の経度
y = 1;
for (i = 0; i < 320; i++) {
x = 1;
for (j = 0; j < 320; j++) {
// 処理対象のピクセル → [X座標=変数x, Y座標=変数y]のピクセル
r = /*処理対象のピクセルの色（R／赤）*/
g = /*処理対象のピクセルの色（G／緑）*/
b = /*処理対象のピクセルの色（B／青）*/
a = /*処理対象のピクセルの色（A／透明度）*/
if (/*ピクセルが透明ではないとき*/) {
if (/*rがおよそ250 かつ gがおよそ230 かつ bがおよそ150*/) int = "4";
else if (/*rがおよそ255 かつ gがおよそ230 かつ bがおよそ0 */) int = "5-";
else if (/*rがおよそ255 かつ gがおよそ153 かつ bがおよそ0 */) int = "5+";
else if (/*rがおよそ255 かつ gがおよそ40 かつ bがおよそ0 */) int = "6-";
else if (/*rがおよそ165 かつ gがおよそ0 かつ bがおよそ33 */) int = "6+";
else if (/*rがおよそ180 かつ gがおよそ0 かつ bがおよそ104*/) int = "7";
North = lat2 + ((lat - lat2) / 320) * i;
South = lat2 + ((lat - lat2) / 320) * (i + 1);
West = lng + ((lng2 - lng) / 320) * j;
East = lng + ((lng2 - lng) / 320) * (j + 1);
/*「メッシュの配列」に追加*/({
type: "Feature",
properties: { rgb: [r, g, b], Intensity: int },
geometry: {
type: "Polygon",
coordinates: [
[
[West, North],
[East, North],
[East, South],
[West, South],
[West, North],
],
],
},
})
}
if (/*変数jが偶数なら*/) x = x + 3;
else x = x + 2;
}
if (/*変数iが偶数なら*/) y = y + 2;
else y = y + 3;
}
geojson = {
type: "FeatureCollection",
features: /*メッシュの配列*/,
}
//hogehoge /*hogehoge*/ ←コメント
かえって分かりにくい説もあり。
ちゃんとできてるの？
生成したのデータが少しでも不正確では、アプリで使用できません。成果物を確認してみましょう。
正しく変換できているのかを確認するには、比較用に正しいデータが必要です（あたりまえ）。しかし、気象庁HPの画像は、不鮮明で確認には向きません。そこで！
こんなサイトがあります。「Jamdp」というサイトです。ご存知の方も多いかもしれませんが、このサイトで、推計震度分布を確認できます。画像ではなくベクターで確認できるため、今回の比較に便利です。これを利用して確認してみましょう。
比較
左：今回の成果物 中央：気象庁HP（今回のデータ元） 右：Jamdp（参考）
３つの方法で、同じデータを表示してみました。現時点では同じに見えますね。では、拡大してみましょう。
気象庁HPのデータは、拡大するとぼやけてしまいました。しかし、その気象庁HPの画像をもとに作成した今回の成果物はGeoJSON形式であるため、ぼやけません。また、今回の成果物を左側の正規のデータと見比べると、誤りなく再現できていることがわかります。また、上の例では、今回の成果物について、震度を独自の配色で表現しています。
まとめると、今回の成果物は、本来のデータを完全に再現できていると言ってよいと思います。
ちなみに、現在は250mメッシュ1で提供される推計震度分布ですが、2023年1月以前は1kmメッシュ1で提供されていました。このプログラムに、1kmメッシュの推計震度分布図の画像を読み込ませても、正常に動作することを確認しました。（メッシュ１個がポリゴン16個で表されるため、非常に非効率ではありますが）
変換したメリット
① 色の変更が容易
推計震度分布図では、震度７は紫、６強は赤...という風に、色と震度階級を対応させます。気象庁ホームページの画像では、気象庁独自の配色が使用されますが、その画像をそのまま自作アプリで使うと、アプリ内の配色と矛盾してしまうことがあります。
そこで、色の変更をしようとするのですが、画像形式では、以下のように、震度階級に属さない「微妙な色」が発生してしまい、変更は困難でした。
しかし、今回の成果物はGeoJSON形式。色の変更はもちろん、パターンで塗りつぶすなど、自由に見た目を変更できます。
② 境界が鮮明
画像形式の推計震度分布図は、ズームしていくと、境界がぼやけてしまいます。見た目があまりよろしくないだけではなく、自分の地域の震度を確認する際、境目が見にくく不便、といったこともあるかもしれません。
しかし、GeoJSON形式はベクターデータであるため、このような問題がなく、鮮明な推計震度分布図が得られます。
推計震度分布図は、揺れの広がりを面的に表現するものですから、境目が見えるか否かにそもそも意味はありませんが...やはり境界がぼやけていると気になりますよね。
③ データ処理
画像データとは違い、GeoJSON形式は数値的なデータを含めることができます。数値データは、様々なことに利用できます。下にその一例を示しておきます。
- 特定の地点の震度を取得する
- データをもとに高度な作図をする
- 統計学的分析をする
- 機械学習の学習材料にする？？
おわりに
今回は、気象庁ホームページの画像から、数値データであるGeoJSONに変換しました。記事を書き始めた当初は、画像がらデータを完璧に復元することができる確証もなかったのですが、結果的に完璧な生成に成功しました。
「無料で推計震度分布図のGeoJSONがほしい」という一見ニッチなニーズは、案外根強くあるようなので、この記事が皆様のお役に立てばうれしいです。
「Jamdp」の掲載を快諾いただいたはち様、ありがとうございました。
宣伝
現在開発中のアプリ、Zero Quakeにも、推計震度分布図の表示機能を搭載しています。
国内外の地震情報、津波情報、緊急地震速報の表示の他、地震検知機能にも対応した、多機能かつシンプルなアプリです。ぜひ、ご利用ください。
また、当アプリはGithubにてオープンソースで公開しております。推計震度分布図の実装など、ご参考にどうぞ。
ちなみに、当アプリでは、パフォーマンスを考慮して、推計震度分布図を GeoJSONではなく、320×320pxのPNG画像に変換して、地図に張り付けています。