場所の緯度経度を取得してマップに統計データを重ねる

はじめに

地理情報を含むデータ分析の際、データを地図上に重ねるとわかりやすいです。
Pythonライブラリのgeocoder、leafmapを用いると簡単に可視化できます。

データ取得

警察庁のWebサイトから、平成30年の「都道府県別刑法犯の認知件数、検挙件数、検挙人員」がまとめられたcsvデータをダウンロードしました。
リンク：https://www.npa.go.jp/hakusyo/r01/data.html

そのうち、九州で起きた犯罪を都道府県別に取得しました。
コードは以下のようになります。
geocoderで都道府県名から緯度と経度を取得しています。

import geocoder
import pandas as pd

df =pd.DataFrame({'都道府県':['福岡県', '佐賀県', '長崎県', '熊本県', '大分県', '宮崎県', '鹿児島県', '沖縄県'],
                '刑法犯総数(認知件数)':[36701, 3581, 3622, 6932, 3331, 4205, 6704, 6878]})


df['緯度経度'] = df['都道府県'].map(lambda x: geocoder.osm(x, timeout=5.0).latlng)
df['緯度'] = df['緯度経度'].map(lambda x: x[0])
df['経度'] = df['緯度経度'].map(lambda x: x[1])

コードを実行すると、データフレームが以下のように取得できます。

可視化コード

import geocoder
import leafmap

# マップの初期表示範囲を指定
location = '九州'
ret = geocoder.osm(location, timeout=5.0)
m = leafmap.Map(center=ret.latlng, zoom=7)

# 犯罪件数に比例した半径の円を追加していく
for i, item in df.iterrows():
    m.add_circle_markers_from_xy(df.iloc[[i]], x='経度', y='緯度', radius=int(item['刑法犯総数(認知件数)']/1000))

m

その他、leafmapでは、航空写真も利用できます。

mapの活用例

GDPやincome情報を含む世界地図

m = leafmap.Map(center=[0, 0], zoom=2)
url = "https://raw.githubusercontent.com/giswqs/leafmap/master/examples/data/countries.geojson"
m.add_vector(
    url, layer_name="Countries", fill_colors=['red', 'yellow', 'green', 'orange']
)
m

地価の可視化

import numpy as np
import matplotlib.colors as cl
import matplotlib.pyplot as plt
import geopandas as gpd
import folium


path_shp = 'L01-22_26_GML/L01-22_26.shp'

gdf = gpd.read_file(path_shp, encoding='shift-JIS')

gdf = gdf[['L01_006', 'geometry']]

# マップの中心を決める（手打ちで丁度よい緯度経度を探しても良い）
y_mass, x_mass = np.array([[geo.y, geo.x] for geo in gdf["geometry"]]).mean(axis=0)

# マップの作成
map1 = folium.Map(location=[y_mass, x_mass], zoom_start=10)


# カラーマップと対応する数値の範囲を設定する
norm = cl.Normalize(vmin=np.log(gdf["L01_006"]).min(), vmax=np.log(gdf["L01_006"]).max())


# カラーマップとして viridis （青→黄）を使用
convert = plt.get_cmap('jet')

# map1 の上に点を打っていく
for geo, pri in zip(gdf["geometry"], np.log(gdf["L01_006"])):
    folium.Circle(location=(geo.y, geo.x),
        radius = 100,
        color = cl.to_hex(convert(norm(pri))),
        fill = True 
    ).add_to(map1)


# map1 を表示
map1

まとめてダウンロード
国土数値情報のファイルを一個一個クリックしてダウンロードするのが面倒な場合は、chromeディベロッパーツールにjavascriptで数行のコードを書くとダウンロードできます。

// 「ダウンロードしますか？」というポップアップが出ないようにする
// https://www.tottio.net/auto-confirm-dialog/
var realConfirm=window.confirm;
window.confirm=function(){
window.confirm=realConfirm;
return true;
};

// xpath を document.evaluateで読み込んで、ループでデータダウンロード
// document.evaluateではなく$x()を使うと、awaitとの併用で $x() is not defined というエラーが出てしまう
for (var i=3; i<50; i++){
console.log(i+'番目');
await new Promise(resolve => setTimeout(resolve, 1000));
document.evaluate('/html/body/main/div/ul['+i+']/li/div[2]/table[2]/tbody/tr[41]/td[6]/a', document).iterateNext().click();
}

全国の地価

バブル期や歴史的な推移も気になります。

慶長５年（１６００）、徳川家康は関ヶ原の戦いに勝つと、天下統一のために全国の街道整備に着手し、翌慶長６年（１６０１）に東海道に宿駅伝馬制度をしきました。
宿駅伝馬制度とは、街道沿いに宿場を設け、公用の旅人や物資の輸送は無料で次の宿駅まで送り継ぐという制度です。輸送のために必要な人馬は、宿場が提供するというものです。
輸送の範囲は原則として隣接する宿場までで、これを越えて運ぶことは禁止されていました。したがって人足と馬もそこで交替することになり、隣の宿場に着くと荷物を新しい馬に積み替えることになります。
東海道には江戸から京都までの間に５３の宿場がありましたから、江戸から京都まで運ぶ場合、５３回の継ぎ替えをすることになります。そのため俗に「５３次」と呼ばれるようになりました。

lineデータの表示

import geocoder
import folium
import geopandas as gpd

# 海岸データの読み込み
path_shp = 'C23-06_01_GML/C23-06_01-g_Coastline.shp'
gdf = gpd.read_file(path_shp, encoding='shift-JIS')


# マップの初期表示範囲を指定
location = '北海道'
ret = geocoder.osm(location, timeout=5.0)
m = folium.Map(location=ret.latlng, zoom_start=8)

# 可視化のためにlineデータ追加
for item in gdf['geometry']:
    points = []
    for lat, lng in zip(item.xy[1], item.xy[0]):
        points.append([lat, lng])
    folium.PolyLine(points).add_to(m)

m

Polygon
foliumを用いてpolygonを表示するにはこのようにします。
leafmapはおそらくfoliumのwrapperな気がします。
そのため、foliumの方が拡張性は高そうです。

import folium
import geopandas as gpd

# 海岸データの読み込み
path_shp = "P21-12_01_GML/P21-12a_01.shp"
df = gpd.read_file(path_shp, encoding='shift-JIS')

m = folium.Map(location=[37, 135], zoom_start=8)


for _, r in df.iterrows():
    # Without simplifying the representation of each borough,
    # the map might not be displayed
    sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
    geo_j = sim_geo.to_json()

    geo_j = folium.GeoJson(data=geo_j,
                           style_function=lambda x: {'fillColor': 'orange'})
    folium.Popup(r['P21A_002']).add_to(geo_j)
    geo_j.add_to(m)
m