ルー○柴が満面の笑みでダイブしてきそうなタイトルですが、やっていきますMIERUNE Advent Calendar 2020 21日目。
境界データの入手
all japanのboundary dataはe-Statからvery easyにgetすることができます。
boundary dataは以下のリンクから色々searchしてgetすることができます。
- 都道府県・市区町村ごと
- shp/XML/GMLの3種類
- 世界測地系緯度経度 or 世界測地系平面直角座標系
ダウンロードしてみる
心の中の大柴がなにやら囁いていましたがスルーしていきましょう。
以下のURLで北海道のshpファイルが入手できます。
クエリパラメーターの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最高ですね！！！！！
皆さんもどんどん利用していきましょう！