地域メッシュとタイルのグリッド比較

## はじめに

最近、地図タイル（ラスタタイル）をメッシュデータとして扱うアイデアを転がしています。

https://qiita.com/mg_kudo/items/62498b053d841644c4e2

日本でメッシュデータと言えば、総務省統計局の「地域メッシュ」があります。統計局が実施する国勢調査の他、さまざまな情報がこの地域メッシュの規格に基づき提供されています（例：国土数値情報、基盤地図情報等）。

https://www.stat.go.jp/data/mesh/m_tuite.html

## 地域メッシュと地図タイルの大きさ比較

地域メッシュは、要は一定の経緯度間隔で区切られたほぼ正方形のグリッドとなります。経緯度で区切っているため、緯度に応じて大きさは異なりますが、おおよその大きさは以下の通りとなります。

なお、[以前の記事](https://qiita.com/mg_kudo/items/9f0b7a4f108ad96d9c25)でまとめたとおり、地図タイルについても１ピクセル当たりのおよその大きさを計算することができますので、地域メッシュの各区画と対応するズームレベル（ZL）も併記してみます。


||地域メッシュの<br>区画の種類|１辺の<br>およその大きさ|対応する<br>タイルの ZL※|
|:--:|:--:|:--:|:--:|
|1|第1次地域区画|約 80 km|ZL 0～1|
|2|第2次地域区画<br>（統合地域メッシュ）|約 10 km|ZL 3～4|
|3|基準地域メッシュ<br>（第3次地域区画）|約 1 km|ZL 7|
|4|2分の1地域メッシュ<br>（分割地域メッシュ）|約 500 m|ZL 8|
|5|4分の1地域メッシュ<br>（分割地域メッシュ）|約 250 m|ZL 9|

※日本周辺の１ピクセルの大きさが各地域メッシュに近いもの（[参考](https://qiita.com/mg_kudo/items/9f0b7a4f108ad96d9c25#%E3%82%BF%E3%82%A4%E3%83%AB%E5%91%A8%E3%82%8A%E3%81%AE%E6%95%B0%E5%AD%97)）

## 地域メッシュをタイルへ変換する

地域メッシュで提供されているデータをタイルデータへ変換するにはどのようにすればよいでしょうか？ポイントは以下の２点あるように思われます。

1. 変換のためのアルゴリズム
2. 各区分の種類に応じた適切なズームレベル

### 変換アルゴリズム

メッシュを変換するのは、結構めんどうくさいです（だからこそ、国が地域メッシュという規格を作り、同じ区画でデータを蓄積し続けることに価値があるのですが）。いくつか変換の際のアルゴリズムを挙げてみたいと思います。

#### 面積按分

メッシュデータを変換する場合は、グリッドの面積に応じて重みづけを行い、按分するのが一般的かと思います。ChatGPT に聞いたところ、[turf.js](https://turfjs.org/) を用いて以下のようなコード（一部修正）を提案してくれました。ただ、こちらの方法はなかなか遅いです。

```JavaScript:面積按分
/**
 * グリッドデータを変換する関数
 * @param {Object} sourceGrid - 元のグリッドデータ（GeoJSON形式）
 * @param {Object} targetGrid - 変換後のグリッドデータ（GeoJSON形式）
 * @param {String} valueProperty - 按分する値のプロパティ名（例: "population"）
 * @returns {Object} - 按分されたターゲットグリッド（GeoJSON形式）
 */
function convertGrid(sourceGrid, targetGrid, valueProperty) {
  // 各ターゲットグリッドのポリゴンに按分された値を計算
  targetGrid.features.forEach((targetFeature) => {
    let totalValue = 0;
    // ターゲットグリッドと交差するソースグリッドを取得
    sourceGrid.features.forEach((sourceFeature) => {
      const polys = turf.featureCollection([targetFeature, sourceFeature]);
      let intersection;
      try {
        intersection = turf.intersect(polys);
      } catch (exceptionVar) {
        //console.log(exceptionVar);
      }
      
      if (intersection) {
        // 交差部分の面積を計算
        const intersectionArea = turf.area(intersection);
        const sourceArea = turf.area(sourceFeature);
        // 面積按分
        const sourceValue = sourceFeature.properties[valueProperty] || 0;
        const apportionedValue = (sourceValue * intersectionArea) / sourceArea;
        // ターゲットグリッドの合計値を更新
        totalValue += apportionedValue;
      }
    });
    // ターゲットグリッドに按分された値を設定
    targetFeature.properties[valueProperty] = totalValue;
  });
  return targetGrid;
}
```

#### 平均

他の方法として、変換先のグリッドに最も近い点を採用する方法もあるかと思います。総数が決まっている人口等については、この方法だと齟齬が出てきますが、例えば、標高等の内挿が有効なデータについては、こちらも採用の余地があります。

例えば、国土地理院の[標高タイル](https://maps.gsi.go.jp/development/hyokochi.html)は、元となる標高モデルデータの標高点のうち、標高タイルの各ピクセルの中心位置の座標に最も近い４つの標高点の値を線形的に平滑化することにより、標高タイルの各ピクセルの中心位置における標高値を算出しています。

#### 値をそのまま引き継ぐ

もっと単純に、グリッドの値をそのまま引き継ぐという方法も、場合によっては使うことになるかもしれません。今回は、変換先のグリッド（タイルの各ピクセル）の代表点が含まれる変換前のグリッドの値を引き継ぐようなコードを作成してみました。

```JavaScript:値をそのまま引継ぐ
function convertGrid(sourceGrid, targetGrid, valueProperty) {
  // 各ターゲットグリッドの値を取得
  targetGrid.features.forEach((targetFeature) => {
    // ターゲットグリッドの代表点が含まれるソースグリッドを取得
    for(let i=0; i<sourceGrid.features.length; i++){
      const sourceFeature = sourceGrid.features[i];
      const point = turf.centroid(targetFeature);
      let intersection;
      try {
        intersection = turf.booleanContains(sourceFeature, point); // return <boolean>
      } catch (exceptionVar) {
        //console.log(exceptionVar);
      }
      
      if (intersection) {
        targetFeature.properties[valueProperty] = sourceFeature.properties[valueProperty];
      }
    }
  });
  return targetGrid;
}
```

### 適切な変換先ズームレベル

地域メッシュの各区分のグリッドから地図タイルへ変換する際は、地域メッシュの大きさに応じて、似ている大きさのピクセルサイズとなるズームレベルを選べば良いかと思われます。

実際、国土地理院の[標高タイル](https://maps.gsi.go.jp/development/hyokochi.html)では、10 m メッシュについては、１ピクセルの大きさが約 8 m の ZL 15 まで、5 m メッシュについては、１ピクセルの大きさが約 4 m の ZL 15 まで作成されています。

一方、地域メッシュは単純に経緯度差が一定のグリッドですが、地域メッシュのピクセルは、Web メルカトルで投影された地図を区切ったものですので、１ピクセルのあたりの経緯度差は緯度が高くなるほど引き延ばされます。ですので、緯度により地域メッシュとタイルピクセルで大きさの比率は異なりますので、特定の地域ばかり見てテストしていると、全国に適用したときに問題が発生する可能性があるため注意です。

|位置|地域メッシュの大きさ<br>（4分の1地域メッシュ）|タイルピクセルの大きさ<br>（ZL 9）|比率|
|:--:|:--:|:--:|:--:|
|那覇<br>（26.15°N 付近）|72121 m²|75007 m²|1.04|
|東京<br>（35.71°N 付近）|65367 m²|61499 m²|0.94|
|札幌<br>（43.04°N 付近）|58941 m²|49917 m²|0.84|

※グリッドの面積は地理院地図の計測機能を用いて算出

また、変換の際は、最大で変換後のグリッドの半分の大きさ分ずれることになります。グリッドのサイズを小さくすればするほど、このずれは小さくなります。ただし、タイルのグリッドを小さくするということは、ズームレベルを深くすることになりますから、タイルの総数がいたずらに増えて[莫大なコスト](https://qiita.com/mg_kudo/items/b3ae4071a69d4026550c)がかかってしまいますので、ある程度のところで妥協する必要があります。

## 実際の比較

地理院地図で提供されている「地域メッシュ統計 人口（総務省統計局）」では、令和2年国勢調査の4分の1地域メッシュ（250mメッシュ）人口が提供されています。地理院地図で使われているデータは GeoJSON ですので、元の地域メッシュの形状は保たれています（ただし、GeoJSON はタイル単位に分割されています）。

こちらを用いて、ためしに 上記の２つのアルゴリズム（面積按分で変換する方法と値をそのまま変換後のメッシュに移す方法）とで比較してみます。後者は変換後のグリッドを集計するとでたらめの総人口となってしまいますので、完全に「見た目」用の変換となります（とはいえ、元データが見た目用のデータなので、一応合理的な割り切りではあります）。

また、変換後のグリッドの大きさが与える影響を見るため、ZL 9 よりもさらに細かい ZL 10～11、緯度の違いを見るため３地域（那覇・東京・札幌）分見てみます。

:::note warn
元になる GeoJSON はタイル境界でグリッドが分割されていますが、含まれている属性値は元となるメッシュ内の人口ですので、タイル境界での分割に対してデータが按分されているわけではありません。もともと表示用のデータだと思いますので、ここでは想定外の利用法をしている点、ご承知おきください。
:::
:::note warn
以下の図の凡例はそれぞれ異なります。QGIS のスタイル機能を用いており、細かい設定は QGIS に丸投げしています。
:::

### 面積按分

* 那覇
![面積按分（那覇）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/a4d38e0c-480f-ec5f-007c-74f254e524f4.png)
* 東京
![面積按分（東京）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/07050871-61c3-4b11-1b13-d07f79520e8d.png)
* 札幌
![面積按分（札幌）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/eadcdf92-8b29-3562-f0c3-0d369ac91049.png)

### 値をそのまま引継ぐ（見た目用）
* 那覇
![値引継ぎ（那覇）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/6d58f0d1-d88d-4d0b-bac2-c6cc3b0e1b56.png)
* 東京
![値引継ぎ（東京）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/d4b82c55-d5f5-0166-03e2-6ba599616912.png)
* 札幌
![値引継ぎ（札幌）.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/3e631cc1-205a-801a-01c2-d409613481e8.png)

## レポジトリ

利用したコードやサンプルデータは以下のレポジトリに置いています。

https://github.com/mghs15/meshToTileGrid

## 感想

面積按分であれば、変換先ズームレベルは、実態に応じた値（＝色）になるので、地域メッシュの大きさに応じて、似ている大きさのピクセルサイズとなるズームレベルで十分な気がします。見た目上は、アンチエイリアスがかかったような感じになりますね。

標高のような平均値を用いる方法については試していませんが、面積按分と同様の考えとなるような気がします。

一方、値をそのまま引き継ぐ場合は、値が同じ（＝色が同じ）であることから、グリッドがずれたような印象になるので、ひとつ分 ZL を深く作っても良いかなぁという感じがします。

面積按分の場合、変換処理にかなり時間がかかってしまうので、面積按分なら ZL を浅くするのは合理的だと思います。ただ、それでも全国分のデータを回すのはしんどそうです。

メッシュ統計とタイルの活用は、今後もいろいろとアイデアを試していきたいと思います。