📚

国土を1ファイルにする — FlatGeobufという選択

に公開

GISや、防災データを扱っていると、「シェープファイルの管理がめんどくさい...」だったり、そもそも「GISのファイルが欠けて読めなくなった」ということがあります。
国土数値情報はデータとしては価値が高い一方で、運用しずらいという課題もあります。
国土数値情報のデータを例に、高速化つシンプルに扱えるようにする方法について解説します。

国土数値情報生データ

国土数値情報ダウンロードサイトとは、国土交通省が公開してる、
日本全国の地理・インフラ・防災データを無料でダウンロードできる公式サイトです。

https://nlftp.mlit.go.jp/ksj/index.html

FlatGeobuf(fgb)

https://flatgeobuf.org/

地理空間データ用に設計されたバイナリフォーマットのことで、GISデータに対してアクセスできる読み込みが非常に速いことで有名です。

通常GISデータは、1つのデータセット(シェープ形式)が、.shp,.dbf,.shx,.prjの4ファイルに分かれていて、4つのファイルのうち一つでも欠けてしまうと正しく読み込めなくなります。fgbはこれらを1つのファイルとして管理できます。

現状では、国土数値情報生データにおいては、FlatGeobuf形式では、提供されていません。
今回は、洪水浸水想定区域(河川単位) 2024年度(令和6年度)版を使用して、fgbを作成していきたいと思います。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31a-2024.html

シェープ形式の場合

河川が6本ある場合、24ファイルになる。
河川A.shp  河川A.dbf  河川A.shx  河川A.prj
河川B.shp  河川B.dbf  河川B.shx  河川B.prj
河川C.shp  河川C.dbf  河川C.shx  河川C.prj
河川D.shp  河川D.dbf  河川D.shx  河川D.prj
河川E.shp  河川E.dbf  河川E.shx  河川E.prj
河川F.shp  河川F.dbf  河川F.shx  河川F.prj
ファイル 内容 欠落時の影響
.shp ジオメトリ座標 データ自体が開けない
.dbf 属性テーブル 図形は表示できるが属性情報がない
.shx 空間インデックス ランダムアクセス不可、読み込み低速化
.prj 座標参照系定義 CRS 不明で重ね合わせ不可

※CRSとは、地図の測り方のルールを表します。日本測量だと、EPSG:6668を使用し、Web/GPSだと、EPSG:4326を使用したりします。(※ EPSG: European Petroleum Survey Group)

fgbの場合

管理が簡単になる。

河川A.fgb
河川B.fgb
河川C.fgb

シェープ形式をfgbに変換

使用するモジュールをインストールします。

pip install pyogrio

新しくディレクトリを作り下記のように配置します。

.
├── data
│   ├── A31a-20-24_24_2400010001_10.dbf
│   ├── A31a-20-24_24_2400010001_10.prj
│   ├── A31a-20-24_24_2400010001_10.shp
│   └── A31a-20-24_24_2400010001_10.shx
└── infer.py

infer.py

import glob
import pyogrio

# .shx が無い場合に .shp から復元する
pyogrio.set_gdal_config_options({"SHAPE_RESTORE_SHX": "YES"})

# data/ 内の shapefile を fgb 化して表示
for shp in glob.glob("data/*.shp"):
    gdf = pyogrio.read_dataframe(shp)
    out = shp.replace(".shp", ".fgb")
    pyogrio.write_dataframe(gdf, out, driver="FlatGeobuf")
    print(f"書き出し: {out}")
    print(gdf)
    print(f"CRS: {gdf.crs}")
    print(f"行数: {len(gdf)}")
    print()

実行結果

書き出し: data/A31a-20-24_24_2400010001_10.fgb
         A31a_201 A31a_202 A31a_203 A31a_204  A31a_205                                           geometry
0      2400010001      員弁川       24      三重県         1  POLYGON ((136.57525 35.06858, 136.57525 35.068...
1      2400010001      員弁川       24      三重県         1  POLYGON ((136.57788 35.06346, 136.57788 35.063...
2      2400010001      員弁川       24      三重県         1  POLYGON ((136.5755 35.06938, 136.5755 35.06942...
3      2400010001      員弁川       24      三重県         1  POLYGON ((136.57538 35.07408, 136.57538 35.074...
4      2400010001      員弁川       24      三重県         1  POLYGON ((136.57531 35.07433, 136.57531 35.074...
...           ...      ...      ...      ...       ...                                                ...
17875  2400010001      員弁川       24      三重県         4  POLYGON ((136.48712 35.20479, 136.48712 35.204...
17876  2400010001      員弁川       24      三重県         4  POLYGON ((136.48775 35.205, 136.48775 35.20504...
17877  2400010001      員弁川       24      三重県         4  POLYGON ((136.48475 35.20917, 136.48475 35.209...
17878  2400010001      員弁川       24      三重県         4  POLYGON ((136.51594 35.1435, 136.51594 35.1435...
17879  2400010001      員弁川       24      三重県         4  POLYGON ((136.52512 35.14563, 136.52512 35.145...

[17880 rows x 6 columns]
CRS: EPSG:6668
行数: 17880

作成されたfgbで読み込んでみる

import pyogrio

FGB_PATH = "data/A31a-20-24_24_2400010001_10.fgb"
gdf = pyogrio.read_dataframe(FGB_PATH)
print(f"=== {FGB_PATH} ===")
print(f"CRS: {gdf.crs}")
print(f"行数: {len(gdf)}")
# Geometry(geopandas では geometry)を表示
print("\n--- geometry ---")
print(gdf.geometry)
# A31a_205 列があれば表示
if "A31a_205" in gdf.columns:
    print("\n--- A31a_205 ---")
    print(gdf["A31a_205"])
else:
    print("\n利用可能な列:", list(gdf.columns))

実行結果

=== data/A31a-20-24_24_2400010001_10.fgb ===
CRS: EPSG:6668
行数: 17880

--- geometry ---
0        POLYGON ((136.70625 35.02317, 136.70625 35.023...
1        POLYGON ((136.70575 35.023, 136.70575 35.02304...
2        POLYGON ((136.70581 35.02304, 136.70581 35.023...
3        POLYGON ((136.70538 35.02288, 136.70538 35.022...
4        POLYGON ((136.705 35.02275, 136.705 35.02279, ...
                               ...                        
17875    POLYGON ((136.58069 35.05775, 136.58069 35.057...
17876    POLYGON ((136.58081 35.05771, 136.58081 35.057...
17877    POLYGON ((136.58156 35.0575, 136.58156 35.0575...
17878    POLYGON ((136.58181 35.0575, 136.58181 35.0575...
17879    POLYGON ((136.58194 35.0575, 136.58194 35.0575...
Name: geometry, Length: 17880, dtype: geometry

--- A31a_205 ---
0        1
1        1
2        1
3        1
4        1
        ..
17875    1
17876    1
17877    1
17878    1
17879    1
Name: A31a_205, Length: 17880, dtype: int32

最後に

株式会社 WHERE (旧: 株式会社 Penetrator) は、シリーズ A ラウンドにおいて総額 5.5 億円の資金調達を実施し、不動産テック業界における更なる成長を目指して採用活動を一層強化しています。

エンジニア, デザイナー, カスタマーサクセス, BizDev, 営業, マーケティングなど、事業拡大を支える多様なポジションで共に挑戦していただける方を待っています!!

▽ 会社のカルチャーを知りたい方はこちら

https://www.wantedly.com/companies/company_9924832

▽ 募集職種を知りたい方はこちら

https://www.wantedly.com/companies/company_9924832/projects

Discussion