anorithのブログ

主にGIS系の記事を書く。

登記所備付地図データをGeoPackageにしてQGISで読み込む

QGIS

概要

  • 登記所備付地図をGeoPackageにしてQGISで読み込む
  • 変換はpython（mojxml2geojson, geopandas）を使用する

参考

DLしたxmlをgeojsonに変換

公式からデータをDLする。

Python環境にmojxml2geojsonをインストールして、以下のスクリプトでDLしたzipを解凍してgeojsonに変換する。 （ブログに載せるほどのスクリプトでもないが・・・）

import os
import glob
import shutil
import subprocess

input_file = R"D:\Geo\touki\47201-3600.zip"
out_dir = R"D:\Geo\touki_out"

# DLしたzipを解凍
input_name = os.path.basename(input_file).split(".")[0]
temp_dir = os.path.join(os.path.dirname(input_file), input_name)
shutil.unpack_archive(input_file, temp_dir)

# 中身のzipを解凍
zip_path_list = glob.glob(os.path.join(temp_dir, "*.zip"))
xml_dir = os.path.join(out_dir, "xml", input_name)
for zip_path in zip_path_list:
    shutil.unpack_archive(zip_path, xml_dir)

# xmlをgeojsonに変換
xml_path_list =glob.glob(os.path.join(xml_dir, "*.xml")) 
for idx, xml_path in enumerate(xml_path_list):
    print(f"[{idx+1}/{len(xml_path_list)}] {os.path.basename(xml_path)}")
    cmd = ["mojxml2geojson", xml_path]
    subprocess.run(cmd)
  • xml -> geojsonの変換は結構遅い。大量にやるなら並列化したいところ。

geopandasでgeojsonをGeoPackageに変換

以下スクリプトでDLしたzipの単位でGeoPackageにまとめて保存する。

import os
import glob
import shutil
import geopandas as gpd
import pandas as pd

input_dir = R"D:\Geo\touki_out\xml\47201-3600" # 先ほどgeojsonに変換したディレクトリを指定
out_dir = R"D:\Geo\touki_out\gpkg"
os.makedirs(out_dir, exist_ok=True)

input_name = os.path.basename(input_dir)
geojson_path_list = glob.glob(os.path.join(input_dir, "*.geojson")) 
gdf_list = []
for idx, geojson_path in enumerate(geojson_path_list):
    print(f"[{idx+1}/{len(geojson_path_list)}] {os.path.basename(geojson_path)}")
    gdf = gpd.read_file(geojson_path)
    gdf = gdf[gdf["座標系"] != "任意座標系"]
    gdf_list.append(gdf.to_crs("epsg:4326"))
gdf = pd.concat(gdf_list)
out_path = os.path.join(out_dir, f"{input_name}.gpkg")
gdf.to_file(out_path, driver='GPKG', layer=input_name)
  • 大したデータ量でもないのでGeoDataFrameの時点でまとめて一つのGeoDataFrameに してしまっている。大量に処理するなら小分けしてGeoPackageへの保存は追記モードでやった方がよさげ。
  • 任意座標系はあきらめる。
  • .to_crs("epsg:4326")をやっておかないと以下のようなエラーが出る場合がある。
ValueError: Cannot determine common CRS for concatenation inputs, got ['JGD2011', 'WGS 84']. Use `to_crs()` to transform geometries to the same CRS before merging.
  • 本当に少ししかやっていないので、なんかほかにも例外出るようなデータはあるかも・・・。

QGISで表示してみる

大字コードで色分けすると以下のような感じ。

QGISのフィルタでNOT "地番" LIKE '%地区外%'として地区外のデータは除外して表示している。xml -> geojsonの段階でmojxml2geojsonの実行オプションとして--excludeとしておけばそもそも地区外のデータをgeojsonに含まれないようにもできる。

フィルタに関するよさげなリンクメモ。