Skip to content

Instantly share code, notes, and snippets.

@keichan34
Last active February 10, 2023 07:47
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save keichan34/d6e8f283bb5810d6f4aa8d941f9a824c to your computer and use it in GitHub Desktop.
Save keichan34/d6e8f283bb5810d6f4aa8d941f9a824c to your computer and use it in GitHub Desktop.
法務省データもろもろ
cp -r moj_data/ ./
mkdir all_zips
find . -name '*.zip' -maxdepth 1 | xargs -P 16 -I '{}' unzip '{}' -d ./all_zips

ここで all_zips に展開されたzipがある 次は任意座標と公共座標系を分ける

find ./all_zips -name '*.zip' | xargs -P 16 zgrep -l '<座標系>任意座標系</座標系>' > ninni_zahyou.txt
mkdir ignore
cat ninni_zahyou.txt | xargs mv -t ./ignore/

分けた

$ ls all_zips/ | wc -l
136101
$ ls ignore/ | wc -l
180142

つまり公共座標系になっているXMLは136101件

次はGeoJSONに変換

下記を convert_in_place.sh として保存

#!/bin/bash -e

unzip -o "$1" -d "$(dirname "$1")"
xml_file="${1%.zip}.xml"
mojxml2geojson "$xml_file"
rm "$xml_file"
geojson_file="${1%.zip}.geojson"
epsg4326_geojson_file="${1%.zip}.epsg4326geojson"
# generated file is in EPSG:6668, we will convert to 4326 for tippecanoe
ogr2ogr -f GeoJSON "$epsg4326_geojson_file" -t_srs "EPSG:4326" "$geojson_file"
ndgeojson_file="${1%.zip}.ndgeojson"
jq -cr '.features | .[]' "$epsg4326_geojson_file" > "$ndgeojson_file"
rm "$geojson_file" "$epsg4326_geojson_file"

並列処理で全部GeoJSONに変換する

find ./all_zips -name '*.zip' -print0 | parallel -0 -j 16 ./convert_in_place.sh

フィルター・マッピング

GeoJSONをフィルタ掛け、propertiesに必要な情報だけ絞って、 GNU Parallel の出力を全て一つのファイルにまとめる。

下記の内容を polygon_filter_script.jq に保存

if (.properties."地番" | test("^[0-9]")) then
{
  type: .type,
  geometry: .geometry,
  properties: (.properties | {
    "市区町村名",
    "市区町村コード",
    "大字コード",
    "丁目コード",
    "小字コード",
    "予備コード",
    "大字名",
    "丁目名",
    "小字名",
    "予備名",
    "地番"
  } )
}
else
empty
end

下記の内容を xml_polygon_generator.sh に保存

#!/bin/bash -e

layer_name=$(basename "$1" .ndgeojson)
ogr2ogr -f geojsonseq /vsistdout/ "$1" -sql "select st_convexhull(st_collect(geometry)) as geometry, \"市区町村コード\",
\"市区町村名\", \"大字コード\", \"丁目コード\", \"小字コード\", \"大字名\", \"丁目名\", \"小字名\" from \"$layer_name\" GROUP BY '1'" -dialect sqlite

下記の内容を point_filter_script.jq に保存

if (.properties."地番" | test("^[0-9]")) then
{
  type: .type,
  geometry: {
    "type": "Point",
    "coordinates": [.properties["代表点経度"], .properties["代表点緯度"]]
  },
  properties: (.properties | {
    "市区町村名",
    "市区町村コード",
    "大字コード",
    "丁目コード",
    "大字名",
    "丁目名"
  } )
}
else
empty
end

ポリゴンのGeoJSONを作る

find ./all_zips -name '*.ndgeojson' -print0 | parallel -0 -j 16 --line-buffer jq -cr -f ./polygon_filter_script.jq '{}' > ./all_polygons.ndgeojson

XMLポリゴンのGeoJSONを作る

find ./all_zips -name '*.ndgeojson' -print0 | parallel -0 -j 16 --line-buffer ./xml_polygon_generator.sh '{}' > ./xml_polygons.ndgeojson

代表点のGeoJSONを作る

find ./all_zips -name '*.ndgeojson' -print0 | parallel -0 -j 16 --line-buffer jq -cr -f ./point_filter_script.jq '{}' > ./all_points.ndgeojson

ベクトルタイル化

mkdir -p $(pwd)/tmp
# ポイントは z0 - z11 まで
tippecanoe \
  -z11 -Z0 \
  -r1 --cluster-distance=10 \
  --read-parallel \
  -f -o all_points.mbtiles \
  -t $(pwd)/tmp \
  -l moj_regmap ./all_points.ndgeojson

# zipポリゴンは z0-13 まで
tippecanoe \
  -z13 -Z0 \
  -pk \
  --drop-smallest-as-needed \
  -pS \
  --detect-shared-borders \
  --read-parallel \
  -f -o xml_polygons.mbtiles \
  -t $(pwd)/tmp \
  -l moj_regmap ./xml_polygons.ndgeojson

# ポリゴンは z14 - z15 まで
tippecanoe \
  -z15 -B15 -Z14 \
  -pk \
  --drop-smallest-as-needed --drop-fraction-as-needed \
  -pS \
  --full-detail=14 \
  --detect-shared-borders \
  --read-parallel \
  -f -o all_polygons.mbtiles \
  -t $(pwd)/tmp \
  -l moj_regmap ./all_polygons.ndgeojson
  
tile-join -pk -f -o moj_regmap.mbtiles ./all_points.mbtiles ./xml_polygons.mbtiles ./all_polygons.mbtiles
  • -pk → 500kb上限を無視する
  • -pS → "Don't simplify lines and polygons at maxzoom"
  • --detect-shared-borders → 共通点をsimplifyしないように(隙間防止)
  • -t ./tmp → デフォルトでは /tmp などを使いますが、もし違う媒体を使ったりしている場合はローカルを指定するといい
95416882 features, 11990840944 bytes of geometry, 221416 bytes of separate metadata, 77839847 bytes of string pool
tile 12/3501/1738 has 209124 features, >200000
Going to try keeping the biggest 71.73% of the features to make it fit
  99.9%  15/28208/13192
real    76m33.786s
user    695m27.907s
sys     9m14.968s

to pmtiles

env TMPDIR=$(pwd)/tmp pmtiles convert moj_regmap.mbtiles moj_regmap.pmtiles
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment