兵庫県高精度3次元地理空間データを利用する

兵庫県が公開している高精度3次元地理空間データを利用する方法について、具体的な手順をまとめました。

## 1. 公開データをダウンロードする

### 1-1. データ公開サイト

まず、下記のサイトにアクセスします。

- [G空間情報センター](https://front.geospatial.jp/) > [兵庫県 50cmメッシュ 数値地形図ポータル（2021～2022年度）](https://www.geospatial.jp/ckan/dataset/2022-hyougo-geo-potal)
    - [兵庫県 50cmメッシュ 標高ラスター/CS立体図（2021～2022年度）](https://www.geospatial.jp/ckan/dataset/2022-hyougo-geo-cs)
    - [兵庫県 50cmメッシュ DSM（2021～2022年度）](https://www.geospatial.jp/ckan/dataset/2022-hyougo-geo-dsm)
    - [兵庫県 50cmメッシュ DEM（2021～2022年度）](https://www.geospatial.jp/ckan/dataset/2022-hyougo-geo-dem)

:::note
ポータルに「点群データダウンロード方法.pdf」が置いてあります
:::

:::note warn
上記とは別に1mメッシュのデータも公開されています。

- [兵庫県_全域数値地形図_ポータル（2010年度～2018年度）](https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-potal)
    - [兵庫県_全域_標高ラスター／CS立体図（2010年度～2018年度）](https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-cs)
    - [兵庫県_全域DSM（2010年度～2018年度）](https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-dsm)
    - [兵庫県_全域DEM（2010年度～2018年度）](https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-dem)

混同しないように注意して下さい。
:::

## 1-2. ダウンロード

下の図のような種類別データの画面で「データ」をクリックします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12310/f3650bf6-7725-4b57-ac65-8838eb57f643.png)

すると、次のようなダウンロード範囲指定画面になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12310/004c3374-3239-4cc5-bfe6-cca4f63ff2ee.png)

ダウンロードしたい位置の升目をクリックして反転表示させると、右上のリストに選択された範囲が追加されます。

そして「ダウンロード」ボタンを押すと、選んだ範囲のデータをダウンロードすることが出来ます。

:::note warn
複数ファイルの連続ダウンロードに対応したブラウザを使用して下さい。そうでないと、最初の一つしかダウンロードされません。

Chrome や Edge なら大丈夫です。
:::

### 1-3. データの形式と内容

ダウンロードしたデータは `zip` 形式の書庫ですが、その内容はデータ種別によって少し異ります。

#### 1-3-1. 標高ラスタ/CS立体図

標高ラスタ/CS立体図では、`zip` 書庫は次のようになっています。
```
[区画名].zip --- [区画名] フォルダ --+-- alt_[区画名]11.tfw
                                   +-- alt_[区画名]11.tif
                                   +-- alt_[区画名]12.tfw
                                   +-- alt_[区画名]12.tif
                                       （中略）
                                   +-- alt_[区画名]43.tfw
                                   +-- alt_[区画名]43.tif
                                   +-- alt_[区画名]44.tfw
                                   +-- alt_[区画名]44.tif
                                   +-- CS_[区画名]11.tif
                                   +-- CS_[区画名]12.tif
                                       （中略）
                                   +-- CS_[区画名]43.tif
                                   +-- CS_[区画名]44.tif
```

[区画名]の後に続く末尾の2桁の数字は、それぞれ上位の区画を `2x2` に分割したもので、次のような入れ子の命名規則になっています。

```
　　　　1 桁目              2 桁目            末尾の2桁
+---------+---------+   +----+----+   +----+----+----+----+
|         |         |   |  1 |  2 |   | 11 | 12 | 21 | 22 |
|    1    |    2    |   +----+----+   +----+----+----+----+
|         |         |   |  3 |  4 |   | 13 | 14 | 23 | 24 |
+---------+---------+   +----+----+   +----+----+----+----+
|         |         |                 | 31 | 32 | 41 | 42 |
+    3    |    4    |                 +----+----+----+----+
|         |         |                 | 33 | 34 | 43 | 44 |
+---------+---------+                 +----+----+----+----+
```

`alt_` で始まるものは標高データの `GeoTIFF` で、`.tif` と `.tfw` で一組のデータです。そして `CS_` で始まるものはCS立体図の `GeoTIFF` です。従って、ファイル数は `16 x 2 + 16 = 48` になります。

個々の `GeoTIFF` ファイルは、2000 x 1500 ピクセルで、東西 1000m x 南北 750m の領域のデータを表わしています。

座標系は `EPSG:2447 JGD2000 / Japan Plane Rectangular CS V` です。

#### 1-3-2. DSM

`DSM` では、zip 書庫は次のようになっています。
```
[区画名].zip --- [区画名] フォルダ --+-- DSM_[区画名]11.txt
                                   +-- DSM_[区画名]12.txt
                                       （中略）
                                   +-- DSM_[区画名]43.txt
                                   +-- DSM_[区画名]44.txt
```

命名規則は、上記の標高ラスター/CS立体図と同じで、ファイル数は 16 です。

個々の `.txt` は、次のような形式のテキスト・ファイルです。

```
63000.250 -116999.750 229.83
63000.750 -116999.750 230.17
63001.250 -116999.750 230.44
63001.750 -116999.750 230.72
```

- 一行に一地点のデータ
- ヘッダ無し
- スペース区切りで、東西位置（`x`）、南北位置（`y`）、高度（`z`）
- 座標系は `EPSG:2447 JGD2000 / Japan Plane Rectangular CS V` で、単位は m
- 2000 x 1500 地点（東西 1000m x 南北 750m）のデータ

#### 1-3-3. DEM

`DEM` では、zip 書庫は次のようになっています。
```
[区画名].zip --- [区画名] フォルダ --+-- DEM_[区画名]11.txt
                                   +-- DEM_[区画名]12.txt
                                       （中略）
                                   +-- DEM_[区画名]43.txt
                                   +-- DEM_[区画名]44.txt
```
:::note
#### DSM と DEM の違い

- DSM ... 数値表層モデル、Digital Surface Model
    - 建物や樹木、橋などの高さを含めた標高を表示するデータ
- DEM ... 数値標高モデル、Digital Elevation Model
    - 建物や樹木などを取り除いた地表面の高さを表示するデータ
:::

:::note warn
#### 標高ラスタ と DEM

ダウンロードできるデータのうち、標高ラスタの `GeoTIFF` と `DEM` のデータは、実質的には同一のデータです。

`DEM` データは `GIS` で扱いやすい `GeoTIFF` に変換してから利用するのが普通ですが、既に変換済みの `GeoTIFF`（標高ラスタ）が存在する訳ですから、手間を掛けて `DEM` データをダウンロードしたり変換したりする必要はありません。
:::

## 2. GeoTIFF への変換

テキストに書かれた数値データである `DSM` と `DEM` は、`QGIS` 等のソフトウェアで扱いやすい `GeoTIFF` 形式のラスタ・データに変換するのが普通です。

### 2-1. QGIS に読み込んで GeoTIFF にエクスポート

もっとも手軽にできる変換方法としては、`DSM` または `DEM` のテキスト・ファイルをドラッグ・アンド・ドロップで `QGIS` に放り込むという手があります。`QGIS` はこの形式の点群データ・ファイルを自動的にラスタ・レイヤとして読み込んでくれますので、それを `GeoTIFF` としてエクスポートすれば良い訳です。

- 参照 : [300万行のテキストを10倍速く地図化する](https://qiita.com/Yfuruchin/items/9fdc130a350570e032c1)

具体的な手順を整理すると、

1. テキスト・ファイルを `QGIS` デスクトップの地図画面にドラッグ・アンド・ドロップで放り込む
2. 「レイヤ DSM_05NFXXXX_05g の CRS を指定して下さい」というメッセージが出るので、`EPSG:2447 JGD2000 / Japan Plane Rectangular CS V` を選ぶ
3. レイヤ・ツリーで当該レイヤを右クリック > エクスポート > 名前を付けて保存で、`GeoTIFF` として保存

と、非常に簡単です。手作業になりますが、変換にかかる時間も数分以内と非常に高速ですので、ファイル数が少ない場合はこれで十分でしょう。

### 2-2. Python スクリプトで一括変換

変換対象のファイル数が 10 を超えると、さすがに手作業は辛くなりますので、Python のスクリプトを使用します。

```python:xyz2tiff.py
import sys
import os
import glob
import pandas as pd # type: ignore
import numpy as np # type: ignore
import rasterio as rio # type: ignore
from rasterio import transform as riotrans # type: ignore

# 兵庫県50cmメッシュの CRS = JGC2000 / Japan Plane Rectangular CS V
CRS = "EPSG:2447"


# メイン
def main(src_dir:str, dst_dir:str):

    # 出力ディレクトリが無ければ作成
    os.makedirs(dst_dir, exist_ok=True)

    i = 0
    for file_in in sorted(glob.glob(os.path.join(src_dir, "*.txt"))):
        base = os.path.basename(file_in)
        name, _ = os.path.splitext(base)
        file_out = os.path.join(dst_dir, f"{name}.tif")
        i = i + 1
        print(f"{i}: ", end="")
        convert(file_in, file_out)
    print(f"Conversion completed. ({i} files)")


# xyz テキスト・ファイルを GeoTIFF に変換する
def convert(file_in: str, file_out: str):

    print(f"Loading {file_in} ... ", end="")
    xyz = load_xyz(file_in)
    print(f"Making matrix ... ", end="")
    arr, west, south, east, north = xyz2matrix(xyz)
    print(f"Writing {file_out} ... ", end="")
    matrix2raster(file_out, arr, west, south, east, north)
    print(f"done")


# xyz テキスト・ファイルを読み出す
def load_xyz(file_in: str) -> pd.DataFrame:
    return pd.read_table(file_in, names="x y z".split(), sep=r"\s+", index_col=False)


# xyz テキスト・ファイルからマトリックスを作成する
def xyz2matrix(xyz: pd.DataFrame) -> (np.ndarray, float, float, float, float): # type: ignore
    mat = xyz.pivot_table(index="y", columns="x", values="z")
    mat.sort_index(axis="index", ascending=False, inplace=True)
    mat.sort_index(axis="columns", ascending=True, inplace=True)
    # Get origin (upper left corner)
    cellsize_x = abs(mat.columns[1] - mat.columns[0])
    cellsize_y = abs(mat.index[1] - mat.index[0])
    south = mat.index.min() - cellsize_y / 2
    north = mat.index.max() + cellsize_y / 2
    west = mat.columns.min() - cellsize_x / 2
    east = mat.columns.max() + cellsize_x / 2

    arr = np.asarray(mat, dtype=np.float32)

    return arr, west, south, east, north


# マトリックスを GeoTIFF に書き出す
def matrix2raster(
    file_out: str, arr: np.ndarray, west: float, south: float, east: float, north: float
):

    transform = riotrans.from_bounds(
        west=west,
        south=south,
        east=east,
        north=north,
        width=arr.shape[1],
        height=arr.shape[0],
    )

    with rio.open(
        file_out,
        "w",
        # driver="GTiff",
        driver="COG",    # Cloud Optimized GeoTIFF
        height=arr.shape[0],
        width=arr.shape[1],
        count=1,
        dtype=str(arr.dtype),
        crs=CRS,
        transform=transform,
        compress="lzw",
        tiled=True,
    ) as raster:
        raster.write(arr, 1)


if __name__ == "__main__":
    if len(sys.argv) < 3:
        sys.stderr.write("Usage: $ python xyz2tiff.py <src_dir> <dest_dir>\n")
        sys.stderr.write("   eg: $ python xzy2tiff.py ./DSM/XYZ ./DSM/TIFF\n")
    else:
        main(sys.argv[1], sys.argv[2])
```

次のように `xyz2tiff.py` に二つの引数を与えて実行します。最初は、`xyz` テキスト・ファイルが置いてあるディレクトリ、次は `GeoTIFF` を生成するディレクトリです。

```bash
python ./xyz2tiff.py ./DSM/XYZ ./DSM/TIFF
```

:::note
`QGIS` デスクトップをインストールした Windows パソコンであれば、`OSGeo4W Shell` を使って上記のスクリプトを実行できます。
:::

:::note warn
`cmd.exe` または `PowerShell` でも実行は可能ですが、下記の環境変数を登録する必要があります。

| 環境変数 | 内容 |
|:----|:----|
| `Path` (追加) | `C:\OSGeo4W\apps\Python312;`<br>`C:\OSGeo4W\apps\Python312\Scripts` |
| `PYTHONHOME` | `C:\OSGeo4W\apps\Python312` |
| `PYTHONPATH` | `C:\OSGeo4W\apps\Python312\Lib;`<br>`C:\OSGeo4W\apps\Python312\DLLs` |

内容は一例です。実際の環境にあわせて設定して下さい。
:::

:::note
`python -V` でバージョンが表示されなかったり、`Python` のインストールを奨めるアプリ・ストアが起動したりする場合は、Microsoft Store の python エイリアスを無効にすると良いです。

1. スタート > 「アプリの設定」
2. 「アプリ実行エイリアス」と検索
3. `python.exe` や `python3.exe` のスイッチを `OFF` にする
:::

### 2-3. Python スクリプトでマージ

下記の `Python` スクリプトで、複数の `GeoTIFF` をマージして一つの `GeoTIFF` を作成することが出来ます。

```python:merge_tiff.py
import os
import sys
import json
import glob
import rasterio as rio # type: ignore
from rasterio.merge import merge as rmerge # type: ignore

# 兵庫県50cmメッシュの CRS = JGC2000 / Japan Plane Rectangular CS V
CRS = "EPSG:2447"


# メイン
def main(config_file: str):

    with open(config_file, "r") as f:
        config = json.load(f)
        src_dir = config["src_dir"]
        dst_dir = config["dst_dir"]

        # 出力ディレクトリが無ければ作成
        os.makedirs(dst_dir, exist_ok=True)

        dst_file_prefix = config["dst_file_prefix"]

        print(f"Src dir = {src_dir} / Dst dir = {dst_dir}")
        
        for job in config["jobs"]:
            keywords = job["keywords"]
            dst_file = f"{dst_file_prefix}{job["dst_file"]}.tif"
            print(f"Keywords = {keywords} / Dst file = {dst_file}")
            merge(src_dir, dst_dir, dst_file, keywords)


# フィルタで指定したファイルをマージして出力
def merge(src_dir: str, dst_dir: str, dst_file: str, keywords):

    all_files = glob.glob(os.path.join(src_dir, "*.tif"))
    in_files = []
    for kw in keywords:
        for fname in all_files:
            if kw in fname:
                in_files.append(fname)

    print(f"Merging rasters ({len(in_files)} files) ...")
    rasters = [rio.open(fn) for fn in in_files]
    arr, transform = rmerge(rasters)
    print(f"Saving {arr.shape[2]} x {arr.shape[1]} raster to {dst_file} ... ", end="")
    with rio.open(
        os.path.join(dst_dir, dst_file),
        "w",
        driver="COG",
        height=arr.shape[1],
        width=arr.shape[2],
        count=arr.shape[0],
        dtype=str(arr.dtype),
        crs=CRS,
        transform=transform,
        bigtiff=True,
        compress="lzw",
        tiled=True,
        num_threads=os.cpu_count(),
    ) as raster:
        raster.write(arr)
    print(f"done.")


if __name__ == "__main__":
    if len(sys.argv) < 2:
        sys.stderr.write("Usage: $ python merge_tiff.py <config.json>\n")
    else:
        main(sys.argv[1])
```

`merge_tiff.py` は次のように構成ファイルを引数として与えて実行します。

```bash
python ./merge_tiff.py ./config.json
```

構成ファイルは、次のような形式の `json` ファイルにして下さい。

```json:config.json
{
  "src_dir": "./DSM/TIFF",
  "dst_dir": "./DSM/TIFF-64",
  "dst_file_prefix": "DSM_",
  "jobs": [
    {
      "dst_file": "00",
      "keywords": ["MF70", "MF71", "MF80", "MF81"]
    },
    {
      "dst_file": "01",
      "keywords": ["MF72", "MF73", "MF82", "MF83"]
    },
    ...
}
```

- `src_dir` ... ソース・ディレクトリ
- `dst_dir` ... デスティネーション・ディレクトリ
- `dst_file_prefix` ... 出力ファイル名のプレフィックス
- `jobs` ... 複数のマージ・ジョブを設定可能
    - `dst_file` ... 出力ファイル名（拡張子無し）
    - `keywords` ... 入力ファイル名のフィルタ

上記のサンプルで最初のジョブは

- 出力ファイル名 `DSM_00.tif`
- 入力ファイル ... `*MF70*.tif`, `*MF71*.tif`, `*MF80*.tif`, `*MF81*.tif` のいづれかに該当するファイル

という構成になります。

:::note
マージされた GeoTIFF は長大なものになり得ますので、"COG" (Cloud Optimized GeoTIFF) ドライバを使い、`bigtiff:True` の設定で書き出しています。
:::

### 2-4. 参照および謝辞

上記の `xyz2tiff.py` と `tiff_merge.py` は、下記の `GitHub Gist` で公開されているスクリプトを参照して作成したものです。

- [Convert .xyz elevation models to GeoTiff](https://gist.github.com/philippkraft/2da0ab4314dd334463fe0e04985bba32)

著者である [Philipp Kraft](https://gist.github.com/philippkraft) に深く感謝の意を表します。

そして、これらのスクリプトについては、原著者 Philipp Kraft に従い、MIT License の下に置きます。

#### Licence (MIT)
Copyright 2025, Nobuo Kihara

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

### 2-5. グリッド

下の図のようなグリッドを作成しておくと便利です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12310/fd54e9ae-7a15-40ae-89b4-71d404f5f475.png)

作り方は、以下のような感じです。

1. x と y について、最小値 と 最大値を調べる
    - マージ前の `GeoTIFF` の左上のものと右下のもの、または、右上のものと左下のものの2つを調べると、全体の領域の `x` と `y` の最大値と最小値が分る
2. ベクタ > 調査ツール > グリッドを作成
    - グリッド・タイプ ... 長方形（ポリゴン）
    - グリッドの範囲 ... x_min, x_max, y_min, y_max で指定
    - 水平方向の間隔 ... 1000.0
    - 垂直方向の間隔 ... 750.0
3. グリッドに "name" 属性を追加
4. グリッドの "name" 属性を編集
    - 属性計算機を使う
    - "row_index", "col_index" を使って "name" を設定する
    - 一回の計算で済ますのは難しいので、複数回に分ける
        - 最初の2桁目まで、3桁目、4桁目、5桁目、6桁目 を別にする
        - `if (条件, concat("name", 'X'), "name")` のような式を使う
            - 条件に合致すれば、`'X'` を追加、合致しなければ元のまま、と言う意味
            - 条件と `'X'` を変更しながら、1桁ずつ追加していく

最後の段階が面倒くさいし、間違いやすいのですが、この `"name"` こそが肝腎の属性なので、とにかく頑張って下さい。ここで頑張っておくと、後々、助かります。