洪水時の浸水深をボクセル風に立体化する

はじめに

水害のリスクを伝える際に、3D にすることでインパクトを与えるというのはよく聞くアイデアです。

国土交通省のハザードマップポータルサイトでは、全国シームレスに洪水浸水想定区域の浸水深を見ることができます。

この浸水深を立体的にしてみれば、浸水リスクが分かりやすいのではと思いましたので、簡易的に立体化させてみることにします。

注意：これはあくまで実験で、実際のリスクを正確に示すものではありません。

データ

実は、ハザードマップポータルサイトのデータは API として提供されております。

その中で、「洪水浸水想定区域（想定最大規模）」のタイルデータ（ラスタタイル）を利用します。

https://disaportaldata.gsi.go.jp/raster/01_flood_l2_shinsuishin_data/{z}/{x}/{y}.png

注意事項は以下の通り。

本データは、洪水浸水想定区域（想定最大規模）_国管理河川と洪水浸水想定区域（想定最大規模）_都道府県管理河川のデータを統合したものです。
都道府県管理河川につきましては、都道府県の許諾を得てタイルデータの配信を行っているため、一部の都道府県のデータ配信のみとなっております。
なお、東京都の公表図面では浸水深0.1m未満の区域は着色されていないため、本サイトにおいても同様の表現としております。
香川県のデータにつきましては、下記利用条件をご確認ください。
利用条件：提供したデータを利用することによって、利用者及び第三者に発生する直接または間接の損失、損害および障害等について、香川県は一切の責任を負わないこと。

いろいろなデータをがっちゃんこさせたように見えます。

立体化

立体化させるには、色々と手法はあるかと思いますが、今回は、マインクラフトのようなボクセル風な表現にしたいと思います。

方法は簡単で、ラスタタイルのピクセルごとに高さのデータ（今回は浸水深）を付与してポリゴンデータにします。その後、ベクトルタイルへ変換し、Mapbox GL JS の fill-extrusion レイヤで立体的に表示させることにします。

WSL1 (Ubuntu 20.04) 上の Node.js v16.13.2 で実行しています。

1. 画像データを取得する

ラスタタイルの画像データ取得・解析するのには、node-fetch と sharp というライブラリを利用しました。特に何も考えずに npm でインストールします。

node-fetch で、画像データを URL から取得してきて、そのデータを sharp で buffer へ変換して返します。

import sharp from 'sharp';
import fetch from 'node-fetch';

const getImageBuffer = (url) => {

  return fetch( url )
    .then( response => {
      return response.arrayBuffer();
    })
    .then( data => {
      const buf = Buffer.from(data);
      return sharp( buf )
        .raw()
        .toBuffer()
    })
    .then( buf => {
      console.log(`GET ${url}`);
      return buf;
    })
    .catch( e => {
      console.log(`ERROR ${url}`);
      console.log(e);
      return null;
    });
}

2. 浸水深を取得する

「洪水浸水想定区域（想定最大規模）」のタイルデータは、以下の通り色分けされています。

画像形式で凡例は出ているのですが、rgb の数値自体は明示されていないので、自分で抽出しました。そのため、この表が完全に正しいという保証はありません。

浸水深 (m)	hex	r	g	b
<0.5	#F7F5A9	247	245	169
0.5-3	#FFD8C0	255	216	192
3-5	#FFB7B7	255	183	183
5-10	#FF9191	255	145	145
10-20	#F285C9	242	133	201
20<	#DC7ADC	220	122	220

ですので、rgb から、対応する浸水深を返すようにすればよいわけです。今回のデータは幅があるのですが、実験なので、代表値を一つだけ返すようにしました。（実際に水害リスクを表現する場合は、もう少し工夫する必要がありそうです。）

const infoFromColor = ( r, g, b ) => {
  if(r == 247, g == 245, b == 169){
    return 0.5;
  }else if(r == 255, g == 216, b == 192){
    return 3;
  }else if(r == 255, g == 183, b == 183){
    return 5;
  }else if(r == 255, g == 145, b == 145){
    return 10;
  }else if(r == 255, g == 133, b == 201){
    return 20;
  }else if(r == 220, g == 122, b == 220){
    return 30;
  }else{
    return -99;
  }
}

3. ピクセルをポリゴンへ変換する

最終的にポリゴンデータに変換する必要があるのですが、そのためには、タイル内のピクセルの座標から、地理座標の経緯度へ変換する必要があります。

Slippy Map Tilenames の説明には、タイル座標（ z/x/y ）を経緯度へ変換するサンプルコードが掲載されていますのでこれを利用します。タイル内のピクセル座標をタイル座標の「小数点以下」として取り扱い、小数点を含めたタイル座標を経緯度へ変換することで、ラスタタイルの各ピクセルに相当したポリゴンを生成できます。

const tile2long = (x,z) => { return (x/Math.pow(2,z)*360-180); }
const tile2lat  = (y,z) => {
  const n=Math.PI-2*Math.PI*y/Math.pow(2,z);
  return (180/Math.PI*Math.atan(0.5*(Math.exp(n)-Math.exp(-n))));
}

取得した画像の各ピクセルごとに色から浸水深を出して、ピクセルの位置からポリゴンを生成し、GeoJSON へ変換し、タイルごとに改行区切りの GeoJSON として出力します。（以下のプログラムでは、直下の dst フォルダへ ${z}-${x}-${y}.ndjson の形式で出力しています。）

import fs from 'fs';

async function img2json ( z, x, y ) {
  const url = `https://disaportaldata.gsi.go.jp/raster/01_flood_l2_shinsuishin_data/${z}/${x}/${y}.png`;
  const buf = await getImageBuffer(url);
  if(!buf) return;
  let res = "";
  const ch = buf.length / ( 256 * 256 ); //今回のタイルは 256 ピクセル四方
  
  // ピクセルごとに演算
  for( let i = 0; i < buf.length / ch; i++ ){
    const [ r, g, b ] = [ buf[ i * ch ], buf[ i * ch + 1 ], buf[ i * ch + 2 ] ];
    
    // RGB を浸水深へ変換
    const dp = infoFromColor( r, g, b );
    
    // ピクセルの位置を経緯度へ変換
    const xx = x + (i % 256) / 256;
    const yy = y + Math.floor(i/256)/256;
    const tx1 = tile2long(xx, z);
    const tx2 = tile2long(xx + 1/256, z);
    const ty1 = tile2lat(yy, z);
    const ty2 = tile2lat(yy + 1/256, z);
    
    // GeoJSON を生成
    let geojson = {
      "type": "Feature",
      "properties": {
        "depth": dp,
        "r": r,
        "g": g,
        "b": b
      },
      "geometry": { 
        "type": "Polygon",
        "coordinates": [[
          [tx1, ty1],
          [tx2, ty1],
          [tx2, ty2],
          [tx1, ty2],
          [tx1, ty1]
        ]]
      }
    };
    
    // 浸水深が 0 m 以上のデータのみ生成
    if(dp > 0){
      res += JSON.stringify(geojson) + "\n";
    }
  }
  
  // ファイルへ書き出し
  fs.writeFile(`./dst/${z}-${x}-${y}.ndjson`, res, (e) => {
    if(e){
      console.log(`ERROR ${z}/${x}/${y} (write file)`);
      console.error(e);
    }
  });
  
  return (`${z}/${x}/${y}`);

}

4. 必要な分だけ 1 から 3 を回す

1 から 3 までの処理が書けたら、必要な分のタイルについて、実際に処理を行います。以下の例は、タイル座標の x が 7270 から 7271、y が 3216 から 3217 までのタイルについて、z=13 のみ取得しています。

あまりに多くのタイルを処理すると、サーバに対して迷惑になるのでご注意を。

最後に、Promise.all() で、処理が終了するかを確認します。

const promiseSet = [];

for( let xi = 7270; xi < 7272; xi++){
  for( let yi = 3216; yi < 3218 ; yi++){
    const pm = img2json(13, xi, yi)
      .catch( e => {
        console.log(`ERROR ${xi} ${yi}`);
        console.log(e);
      });
    promiseSet.push(pm);
  }
}

Promise.all(promiseSet)
  .then( values => {
    console.log(values);
    console.log(`COMPLETED`);
  })
  .catch( e => {
    console.log(`ERROR`);
    console.log(e);
  });

以上の 1 から 4 の処理を一つの js ファイルに記述し、Node.js 上で実行すれば、GeoJSON が生成されていくはずです。

5. タイル化と表示

4 までの処理で、GoeJSON 一式を作成することができたら、tippecanoe を用いてベクトルタイルにします。

コマンドの例は以下の通り。Windows では動かないので、WSL 上にインストールして実行しています。

なお、今回は隣り合うタイルで属性値が同じ場合、融合させるように --coalesce、--reorder、--no-simplification-of-shared-nodes を指定しています。処理が重くなるほか、思った通りに融合してくれない（おそらく、データ側に起因するような気がしているが）ので、ちゃんとやるなら、tippecanoe 任せにせずに、最初に前処理で融合処理をしてしまう方がいいです。

tippecanoe -l suibou -e ./tile/ ./dst/* \
  --force --coalesce --reorder \
  --hilbert --no-simplification-of-shared-nodes \
  --no-tile-size-limit --no-feature-limit\
  --minimum-zoom=12 --maximum-zoom=12 --base-zoom=12 \
  --no-tile-compression --no-line-simplification

上記コマンドでは、dst 内の GeoJSON を変換して、tile というフォルダに出力します。
このタイルを Mapbox GL JS 等で、以下のようにソース、スタイルレイヤを追加して、地図上に表示させれば、ボクセル風な立体化ができます。

map.addSource("shinsui", {
  "type": "vector",
  "tiles":[
    "http://localhost/tile/{z}/{x}/{y}.pbf" // タイルへのパス
  ],
  "maxzoom": 12,
  "minzoom": 12,
  "attribution": "<a href='https://disaportal.gsi.go.jp/hazardmap/copyright/opendata.html#l2shinsuishin' target='_blank'>ハザードマップポータルのデータ</a>を加工"
});
map.addLayer({
  "id": "flood-depth-L2",
  "type": "fill-extrusion",
  "source": "shinsui",
  "source-layer": "suibou",
  "maxzoom": 22,
  "minzoom": 12,
  "filter": [">", ["get", "depth"], 0],
  "layout": {},
  "paint": {
    "fill-extrusion-color": [
      "rgb",
      ["get", "r"],
      ["get", "g"],
      ["get", "b"]
    ],
    "fill-extrusion-opacity": 0.8,
    "fill-extrusion-height": ["get", "alti"]
  }
});

完成品

デモサイト

建物の高さは実際の高さを示していないことに注意

レポジトリ

※上記で紹介したコードは記事用に編集しているので、実際のレポジトリのコードとは違う書き方になっているところがあります。

おまけ

上記コードを応用して、地理院タイルの標準地図を、標高タイルを用いてボクセル風にしてみました。あまり見栄えはよくないかもしれません。最初、標高タイルから読み取った標高値をそのまま表示（小数点以下2桁）していましたが、そうするとなめらかな反面、地理院地図3D の劣化版という印象が強かったため、標高値は 10 m 単位で丸めて、ボクセル感を強調しています。

終わりに

正確な建物データや、もっと細かく浸水深を分類した浸水想定区域図データがあれば、もっと直感的にリスクがわかるような地図ができるのではないかと思います。

リアルな 3D データの方がリアリティがありますが、シンプルでサクサク使えるようなサイトで、気軽に浸水リスクを実感できるというのも、需要があるのではないかと思います。

最後になりますが、本記事は実験ですので、記事や参照レポジトリ等の注意事項を理解いただき、実際のリスク評価に使わないようにしてください。