MLS新規登録数を市町村別に出力

!wget https://nlftp.mlit.go.jp/ksj/gml/data/N03/N03-2021/N03-20210101_38_GML.zip

!pip install geopandas
!pip install rtree
!pip install pygeos

import datetime

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd

url = "https://mls.js2hgw.com/update.php?base_enb=737280"

# 現在の日時
JST = datetime.timezone(datetime.timedelta(hours=+9))
dt_now = datetime.datetime.now(JST).replace(tzinfo=None)

# 日付範囲（８日前から今日まで）
dt_start = (dt_now - datetime.timedelta(days=8)).date()
dt_range = pd.date_range(start=dt_start, end=dt_now.date())

df0 = pd.read_json(url + "&output=json", convert_dates=["created", "updated", "modified"])

# 空白は不明
df0["area1"].mask(df0["area1"] == "", "不明", inplace=True)

# 経度・緯度
df0[["lat", "lng"]] = (
    df0["location"].str.strip("()").str.split(",", expand=True).astype(float)
)

df0["場所"] = df0["pref"].str.cat(df0["area1"])
df0["eNB-LCID"] = df0["enb"].astype(str).str.cat(df0["lcid"].astype(str), "-")
df0["date"] = df0["created"].dt.date

df1 = pd.crosstab(df0["date"], df0["area1"])

df2 = df1.reindex(dt_range, fill_value=0)

# 前日
se = df2.iloc[-2]

if se.sum():

    twit = f"MLS 新規登録\n\n{se.name:%Y/%m/%d}\n\n"

    for k, v in se[se > 0].to_dict().items():
        twit += f"{k} : {v}\n"

    twit += f"\n{url}"

    print(twit)

    print("-" * 30)

    df3 = df0[df0["date"] == se.name].copy().sort_values(["enb", "lcid"])

    # 市町村別eNB-LCID

    cities = df3.groupby(by="area1")["eNB-LCID"].apply(list).to_dict()

    for k, v in cities.items():
        print(f"【{k}】")
        for i in v:
            print(i)
        print()

    geo_df = gpd.GeoDataFrame(df3, geometry=gpd.points_from_xy(df3.lng, df3.lat), crs=6668)
    ehime = gpd.read_file("N03-20210101_38_GML.zip!N03-20210101_38_GML")

    # 位置表示

    base = ehime.plot(color="white", edgecolor="black")
    ax = geo_df.plot(ax=base, marker="o", color="red", markersize=5)
    ax.set_axis_off()

    # グラフを保存
    plt.savefig("map.png", dpi=200)

    plt.show()