Goldilocks Days: 9月 2024

2024年9月30日月曜日

Pythonで地図空間データを扱う③

GeoPandasで都市データのプロット

GeoPandasのdatasetsには、naturalearth_lowresの他に主要都市のデータnaturalearth_citiesがあります。

これを使い世界地図に都市もプロットしてみましょう。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import geopandas as gpd
import matplotlib.pyplot as plt

# 世界地図のデータ
world = gpd.datasets.get_path('naturalearth_lowres')
# 主要都市のデータ
cities = gpd.datasets.get_path('naturalearth_cities')

wdf = gpd.read_file(world)
cdf = gpd.read_file(cities)

print(wdf.head())
print()
print(wdf.shape)
print()
print(cdf.head())
print()
print(cdf.shape)

結果

pop_est continent name iso_a3 gdp_md_est geometry

0 889953.0 Oceania Fiji FJI 5496 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...

1 58005463.0 Africa Tanzania TZA 63177 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...

2 603253.0 Africa W. Sahara ESH 907 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...

3 37589262.0 North America Canada CAN 1736425 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...

4 328239523.0 North America United States of America USA 21433226 MULTIPOLYGON (((-122.84000 49.00000, -120.0000...

(177, 6)

name geometry

0 Vatican City POINT (12.45339 41.90328)

1 San Marino POINT (12.44177 43.93610)

2 Vaduz POINT (9.51667 47.13372)

3 Lobamba POINT (31.20000 -26.46667)

4 Luxembourg POINT (6.13000 49.61166)

(243, 2)

naturalearth_citiesは主要都市の名前と位置というシンプルなデータです。

これを世界地図に乗せていきます。

世界地図のGeoDataFrameオブジェクトを変数にして、axに指定することで両方プロットすることが出来ます。

# ベースとなる世界地図
bg = wdf.plot(column='continent', figsize=(18, 8))
# 主要都市をマッピング
cdf.plot(ax=bg, marker='o', color='firebrick')

# アノテーションで主要都市の名前をマッピング
for i, (city, point) in cdf.iterrows():
    # pointオブジェクトから緯度経度を抽出
    x1, y1, x2, y2 = point.bounds
    bg.annotate(city, xy=(x1, y1))

plt.show()

都市名はannotateでプロットします。

geomatryデータからshapely.geometry.point.Pointオブジェクトを抜き出し、 bounds変数には(minx, miny, maxx, maxy)が入っているので、xy座標に指定します。

shapely.Point

https://shapely.readthedocs.io/en/stable/reference/shapely.Point.html

世界地図に都市名も表示出来ました。

Pythonで地図空間データを扱う②

GeoPandasでコロプレスマップを表示

基本的な動作を確認したら、コロプレスマップを作ってみます。

naturalearth_lowresには大陸以外に人口とGDPデータがあるので、それらを色で段階的に分けて表示します。

人口マップ

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import geopandas as gpd
import matplotlib.pyplot as plt
#日本語を表示するため
import japanize_matplotlib


fn = gpd.datasets.get_path('naturalearth_lowres')
gdf = gpd.read_file(fn)

# 人口データをint型にします
gdf['population'] = gdf['pop_est'].astype("int")
# 人口の多い順にソートします
gdf = gdf.sort_values("pop_est", ascending=False)
gdf = gdf.drop("pop_est", axis=1)

print(gdf.head(10))

continent name iso_a3 gdp_md_est population

139 Asia China CHN 14342903 1397715000

98 Asia India IND 2868929 1366417754

4 North America United States of America USA 21433226 328239523

8 Asia Indonesia IDN 1119190 270625568

102 Asia Pakistan PAK 278221 216565318

29 South America Brazil BRA 1839758 211049527

56 Africa Nigeria NGA 448120 200963599

99 Asia Bangladesh BGD 302571 163046161

18 Europe Russia RUS 1699876 144373535

27 North America Mexico MEX 1268870 127575529

人口を分かりやすく表示しています。

それらのデータをプロットすることで人口のコロプレスマップが出来ます。

gdf.plot(column='population', cmap="jet", legend=True, figsize=(16, 8))
plt.title("人口マップ")
plt.show()

中国、インドの人口が多いのがひと目で理解できます。

各国のGDPマップ

次はGDPデータを使いコロプレスマップにしてみます。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import geopandas as gpd
import matplotlib.pyplot as plt
import japanize_matplotlib


fn = gpd.datasets.get_path('naturalearth_lowres')
gdf = gpd.read_file(fn)

gdf = gdf.rename(columns={"gdp_md_est": "gdp_md"})
# GDPの多い順に並べます
gdf = gdf.sort_values("gdp_md", ascending=False)

print(gdf.head(10))

gdf.plot(column='gdp_md', cmap="jet", legend=True, figsize=(16, 8))
plt.title("GDP")
plt.show()

各国のGDPのコロプレスマップを表示出来ました。

2024年9月27日金曜日

Pythonで地図空間データを扱う

GeoPandasを使いPythonで地図を描画

Pythonはデータ処理関係のライブラリが豊富ですが、地図とデータを組み合わせられたら地域データの詳細もより分かりやすくなります。ここ数年で注目度が高まっているGeoPandasを使って地図データを扱ってみたいと思います。

pip install geopandas

他のライブラリにも依存してしますので適宜インストールしてください。

GeoDataFrameに変換するにはshapeファイルやGeoJSONファイルが必要ですが、テスト用のファイルがdatasetsにあります。

現在は3つほどと少ないですが、気軽にGeoDataFrmeを試すにはとても良いデータです。

naturalearth_lowresを使用してみます。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import geopandas as gpd
import matplotlib.pyplot as plt

# ①地図データファイルのあるパスを指定
path = gpd.datasets.get_path('naturalearth_lowres')

# ②ファイルからGeoDataFrameオブジェクトを生成
gdf = gpd.read_file(path)

print(gdf.head())
print("")
print(gdf.shape)
print("")
print(gdf.crs)

結果

pop_est continent name iso_a3 gdp_md_est geometry

0 889953.0 Oceania Fiji FJI 5496 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...

1 58005463.0 Africa Tanzania TZA 63177 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...

2 603253.0 Africa W. Sahara ESH 907 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...

3 37589262.0 North America Canada CAN 1736425 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...

4 328239523.0 North America United States of America USA 21433226 MULTIPOLYGON (((-122.84000 49.00000, -120.0000...

(177, 6)

EPSG:4326

pandasに慣れていれば非常に扱いやすいですね。

headやshapeも同じ使い方です。

Geoデータフレームはデータにgeometryという位置データが追加されているのがpandasデータフレームとの違いです。

CRSとは座標参照系(Coordinate Reference System)のことで、位置情報のルールを表します。

文字コードのようなものです。

このデータを使ってプロットします。

gdf.plot()
plt.show()

世界地図がプロットされました。 pandasライクで簡単に地図データを扱うことが出来ます。このままでは味気ないので、オプションをつけて、少し装飾してみます。

columnは色分け、cmapはカラーマップを指定します。

gdf.plot(column="continent", cmap="Accent", legend=True, figsize=(16, 8))
plt.show()

綺麗になりました。

2024年9月2日月曜日

2点間を通る直線の式を連立方程式で解く(Python)

2点間を通る直線の式

Pythonで連立方程式を解くことが出来たら、2点間を結ぶ一次方程式も導き出せるので

Pythonで書いてみます。

2点の位置 (-1, -5) (3, 3)

#! /usr/bin/env python  
  
# -*- coding:utf-8 -*-  
   
# 
import numpy as np  
import matplotlib.pyplot as plt  
from sympy import Symbol, solve
  
# 一次関数のデータを作る関数   
def linear(a, b, ran=np.arange(-10,10)):  
  
    arr = []  
    for x in ran:  
        y = x*a + b  
        arr.append(y)  
  
    # 式から配列データを返します
    return arr  
  
ax = fig.add_subplot(1, 1, 1)  
   
# これで正方形になります 
ax.set_aspect('equal', adjustable='box')  
  
# 表示する幅  
ax.set_xlim(-10, 11)  
ax.set_ylim(-10, 11)  
  
# グラフの刻み  
ax.set_xticks(np.arange(-10, 10, 5))  
ax.set_yticks(np.arange(-10, 10, 5))  
   
# 補助線をonにします  
ax.minorticks_on()   
  
# グリッドの主線と補助線  
ax.grid(which="major", color="black", alpha=0.5)  
ax.grid(which="minor", color="gray", linestyle="--", alpha=0.5)  
  
# 数値表示を0に合わせます
ax.spines['bottom'].set_position("zero")  
ax.spines['left'].set_position("zero")  

# xの生成
x = np.arange(-10, 10)

# 2点の位置 (-1, -5) (3, 3)
x1, y1 = (-1, -5)
x2, y2 = (3, 3)

a = Symbol("a")
b = Symbol("b")

# 2点を連立方程式に変えます
s1 = a * x1 + b - y1
s2 = a * x2 + b - y2

result = solve((s1, s2))

# 得られた解
print(result)

# yの生成
y = linear(result[a], result[b], x)

# 得られた一次方程式のプロット
ax.plot(x, y)

# 点(-1, -5) のプロット
ax.plot(x1, y1, marker="o")
ax.annotate("(-1, -1)", xy=(x1, y1), color="red")
# 点(3, 3) のプロット 
ax.plot(x2, y2, marker="o")
ax.annotate("(3, 3)", xy=(x2, y2), color="red")
plt.show()

このように2点間から一次関数のグラフが描画出来ました。

statsmodelsのseasonal_decomposeでエラーが出る場合

seasonal_decompose

statsmodels.apiで波形分解する場合、sm.tsa.seasonal_decomposeを使いますが、バージョンによって書き方が異なります。

TypeError: seasonal_decompose() got an unexpected keyword argument 'freq'

というエラーが出る場合は、新しい書き方に置き換えればエラーは出なくなると思います。

statsmodels==0.14.2

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# sm.tsa.seasonal_decompose(data, freq=365).plot()

# 上記でエラーが出る場合、こちらに変更
sm.tsa.seasonal_decompose(data, period=365, extrapolate_trend='freq').plot()

登録: 投稿 (Atom)