LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@softark

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

Last updated at Posted at 2025-04-05

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

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

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

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

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

1-2. ダウンロード

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

image.png

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

image.png

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

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

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

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

DSM と DEM の違い

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

標高ラスタ と DEM

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

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

2. GeoTIFF への変換

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

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

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

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

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

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

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

変換対象のファイル数が 10 を超えると、さすがに手作業は辛くなりますので、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",
        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 を生成するディレクトリです。

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

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

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

環境変数 内容
Path (追加) C:\OSGeo4W\apps\Python312;
C:\OSGeo4W\apps\Python312\Scripts
PYTHONHOME C:\OSGeo4W\apps\Python312
PYTHONPATH C:\OSGeo4W\apps\Python312\Lib;
C:\OSGeo4W\apps\Python312\DLLs

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

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

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

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

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

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 は次のように構成ファイルを引数として与えて実行します。

python ./merge_tiff.py ./config.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 のいづれかに該当するファイル

という構成になります。

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

2-4. 参照および謝辞

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

著者である Philipp Kraft に深く感謝の意を表します。

そして、これらのスクリプトについては、原著者 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

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

  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" こそが肝腎の属性なので、とにかく頑張って下さい。ここで頑張っておくと、後々、助かります。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?