ぼく「e-StatのboundaryデータをGeoPandasでdissolveしてGeoPackageに保存」

ルー○柴が満面の笑みでダイブしてきそうなタイトルですが、やっていきますMIERUNE Advent Calendar 2020 21日目。

境界データの入手

all japanのboundary dataはe-Statからvery easyにgetすることができます。

boundary dataは以下のリンクから色々searchしてgetすることができます。

https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015&coordsys=1&format=shape

• 都道府県・市区町村ごと
• shp/XML/GMLの3種類
• 世界測地系緯度経度 or 世界測地系平面直角座標系

ダウンロードしてみる

心の中の大柴がなにやら囁いていましたがスルーしていきましょう。

以下のURLで北海道のshpファイルが入手できます。

https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212015&code=01&coordSys=1&format=shape&downloadType=5

クエリパラメーターのcode=01を変更することによってダウンロードする都道府県を変更できます。

都道府県コードはこちらを参照。

スクリプトを書いてみる

ひとまず全都道府県の境界データ（shp）をダウンロードするスクリプトを書いてみましょう。

利用する外部パッケージは以下の二つなのでダウンロードしていきましょう。

• tqdm：プログレスバーを表示してくれる粋なやつ
• requests：簡単にhttpでリクエストできる頼もしいやつ

ちなみに、以下のスクリプトは色々あって対象ファイルの総容量がキャッチできなかったのでtqdmはバーの表示をしてくれませんが、ダウンロード速度は表示してくれます。

pip install tqdm
pip install requests

スクリプトはこんな感じになりました。

import os
import re
import sys
from concurrent.futures.thread import ThreadPoolExecutor
from pathlib import Path

import requests
from tqdm import tqdm

URL = "https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212015&code={pref_code}&coordSys=1&format=shape&downloadType=5"

str_pref = [str(pref_code).zfill(2) for pref_code in range(1, 48)]
args = [(URL.format(pref_code=p), str(Path("./").resolve())) for p in str_pref]


def main():
    """処理を実行します

    """
    with ThreadPoolExecutor() as executor:
        executor.map(lambda p: file_download(*p), args)


def get_file_name_from_response(url, response):
    """responseのContent-Dispositionからファイル名を取得、できなければURLの末尾をファイル名として返す

    Args:
        url (str): リクエストのURL
        response (Response): responseオブジェクト

    Returns:
        str: ファイル名を返す

    """
    disposition = response.headers["Content-Disposition"]
    try:
        file_name = re.findall(r"filename.+''(.+)", disposition)[0]
    except IndexError:
        print("ファイル名が取得できませんでした")
        file_name = os.path.basename(url)
    return file_name


def file_download(url, dir_path, overwrite=True):
    """URLと保存先ディレクトリを指定してファイルをダウンロード

    Args:
        url (str): ダウンロードリンク
        dir_path (str): 保存するディレクトリのパス文字列
        overwrite (bool): ファイル上書きオプション。Trueなら上書き

    Returns:
        Path: ダウンロードファイルのパスオブジェクト

    Notes:
        すでにファイルが存在していて、overwrite=Falseなら何もせず
        ファイルパスを返す

    """
    res = requests.get(url, stream=True)

    parent_dir = Path(dir_path).parent
    file_name = get_file_name_from_response(url, res)
    download_path = parent_dir / file_name

    if download_path.exists() and not overwrite:
        print("ファイルがすでに存在し、overwrite=Falseなのでダウンロードを中止します。")
        return download_path

    # content-lengthは必ず存在するわけでは無いためチェック
    try:
        file_size = int(res.headers['content-length'])
    except KeyError:
        file_size = None
    progress_bar = tqdm(total=file_size, unit="B", unit_scale=True)

    if res.status_code == 200:
        print(f"{url=}, {res.status_code=}")
        print(f"{file_name}のダウンロードを開始します")
        with download_path.open('wb') as file:
            for chunk in res.iter_content(chunk_size=1024):
                file.write(chunk)
                progress_bar.update(len(chunk))
            progress_bar.close()
        return download_path
    else:
        print(f"{url=}, {res.status_code=}")
        print("正常にリクエストできませんでした。システムを終了します。")
        sys.exit(1)


if __name__ == '__main__':
    main()

スクリプトを実行すると13行目、args = [(URL.format(pref_code=p), str(Path("./").resolve())) for p in str_pref] のPath("./")で指定したディレクトリにshpファイルが格納されたzipファイルがダウンロードされます。

ThreadPoolExecutorを利用してシングルスレッドで並行処理を行なっているためダウンロードが爆速です。

案外さらっとall japanのboundary dataをgetすることができましたね！s

GeoPandasで読み込んで結合

ダウンロードしてきたzipファイルをGeoPandasに読み込ませていきましょう。

import glob
from pathlib import Path

import geopandas as gpd
import pandas as pd

zip_path_list = [Path(l) for l in glob.glob("./*.zip")]

town_gdf = pd.concat([
    gpd.read_file("zip://" + str(zipfile))
    for zipfile in zip_path_list
]).pipe(gpd.GeoDataFrame)

town_gdf["AREA_CODE"] = town_gdf['PREF'].str.cat(town_gdf['CITY'])
city_gdf = town_gdf.dissolve(by="AREA_CODE", as_index=False)
pref_gdf = city_gdf.dissolve(by="PREF", as_index=False)

globでzipファイルのリストを作成
→リスト内包表記でgeopandasのGeoDataFrameクラスのリストをzipファイルのリストから作成
→pandas.concatメソッドで縦方向にデータフレームを連結
→最後にgeopandas.GeoDataFrameに変換しています。

今回ダウンロードしたshpファイルは町丁目まで含んだ詳細なポリゴンデータですので、カラムを指定してGeoDataFrame.dissolveメソッドで利用市町村ごと、都道府県ごとのポリゴンも作成しました。

GeoPackageに保存する

はいこれだけ！

pref_gdf.to_file("boundary.gpkg", layer='pref', driver="GPKG")
city_gdf.to_file("boundary.gpkg", layer='city', driver="GPKG")
town_gdf.to_file("boundary.gpkg", layer='town', driver="GPKG")

これでboundary.gpkgというファイル名でpref・city・townの3レイヤ
を持つGeoPackageが作成されました！

表示させてみる

いい感じ！！！！

とても簡単に全国分の境界データを利用できるe-Stat最高ですね！！！！！

皆さんもどんどん利用していきましょう！