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

ベースの地図が出来た所で、他のデータを被せてみます。

国土地理院の 500mメッシュ別将来推計人口データ を使用します。

同じく神奈川県のデータ 500m_mesh_suikei_2018_shape_14.zip をダウンロードします。

ベースの地図データと同じ場所に展開します。

そして同じようにread_fileで読み込みます。

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. #  
5. import geopandas as gpd  
6. import matplotlib.pyplot as plt  
7. import japanize_matplotlib  
8.   
9. fn = "../data/N03-20240101_14.shp"  
10. gdf = gpd.read_file(fn)  
11. fn2 = "../data/500m_mesh_2018_14.shp"  
12. pdf = gpd.read_file(fn2)  
13.   
14. print(gdf.head(10))  
15. print("")  
16. print(gdf.shape)  
17. print("")  
18. print(pdf.head())  
19. print("")  
20. print(pdf.shape)  

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

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

fn = "../data/N03-20240101_14.shp"
gdf = gpd.read_file(fn)
fn2 = "../data/500m_mesh_2018_14.shp"
pdf = gpd.read_file(fn2)

print(gdf.head(10))
print("")
print(gdf.shape)
print("")
print(pdf.head())
print("")
print(pdf.shape)

(6193, 235)と巨大なデータフレームなので、必要なデータ以外は取り除いてもいいでしょう。

view plainprint?

1. # 必要な列だけ取り出し軽量化します。  
2. pdf2 = pdf[['SHICODE', 'PTN_2020', 'geometry']]  
3.   
4. # ベースとなる地図  
5. bg = gdf.boundary.plot(linewidth=0.5, edgecolor="gray", figsize=(14, 8))  
6.   
7. # 人口データのプロット  
8. pdf2.plot(ax=bg, column="PTN_2020", cmap="jet", legend=True)  
9.   
10. # 地域名を表示  
11. for name, row in citys:  
12.     if row['N03_001'].count() > 1:  
13.         x = row["geometry"].centroid.bounds.minx.mean()  
14.         y = row["geometry"].centroid.bounds.miny.mean()  
15.         bg.annotate(name, xy=(x, y), color="red")  
16.     else:  
17.         x = row["geometry"].centroid.bounds.minx  
18.         y = row["geometry"].centroid.bounds.miny  
19.         bg.annotate(name, xy=(x, y), color="red")    
20.   
21. plt.title("神奈川県人口マップ")  
22. plt.show()  

# 必要な列だけ取り出し軽量化します。
pdf2 = pdf[['SHICODE', 'PTN_2020', 'geometry']]

# ベースとなる地図
bg = gdf.boundary.plot(linewidth=0.5, edgecolor="gray", figsize=(14, 8))

# 人口データのプロット
pdf2.plot(ax=bg, column="PTN_2020", cmap="jet", legend=True)

# 地域名を表示
for name, row in citys:
    if row['N03_001'].count() > 1:
        x = row["geometry"].centroid.bounds.minx.mean()
        y = row["geometry"].centroid.bounds.miny.mean()
        bg.annotate(name, xy=(x, y), color="red")
    else:
        x = row["geometry"].centroid.bounds.minx
        y = row["geometry"].centroid.bounds.miny
        bg.annotate(name, xy=(x, y), color="red")  

plt.title("神奈川県人口マップ")
plt.show()

地図データに人口データを重ね合わせることが出来ました。

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

日本地図データ

テストデータで概要を掴めたら、外部からshapeファイルを読み込み使用してみます。

 日本に住んでる方は、日本の地域データを活用することが多いでしょうから国土地理院のデータを読み込んでみます。

国土地理院

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-2024.html

全国はデータが大きいので今回は神奈川県のデータをダウンロードしました。

データはzipファイルなので、使用する場所に移動して解凍します。

$ unzip N03-20240101_14_GML.zip 

解凍するとshapeファイルなどいくつかのファイルが出来ます。

それらをリンクして使用されるので、個別で移動させたり、ファイル名を変えたりはしないでください。

回答したshapeファイルを前回と同じように読み込みます。

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. import geopandas as gpd  
5. import matplotlib.pyplot as plt  
6. import japanize_matplotlib  
7.   
8. fn = "../data/N03-20240101_14.shp"  
9. gdf = gpd.read_file(fn)  
10.   
11. print(gdf.head())  
12. print("")  
13. print(gdf.shape)  
14. print("")  

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

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

fn = "../data/N03-20240101_14.shp"
gdf = gpd.read_file(fn)

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

結果

	N03_001	N03_002	N03_003	N03_004	N03_005	N03_007	geometry
0	神奈川県	None	None	横浜市	鶴見区	14101	POLYGON ((139.65294 35.50000, 139.65281 35.50010, 139.65270 35.50015, 139.65221 35.50030, 139.65214 35.50034, 139.65207 35.50038,・・・
1	神奈川県	None	None	横浜市	鶴見区	14101	POLYGON ((139.67394 35.46229, 139.67347 35.46329, 139.67303 35.46426, 139.67272 35.46494, 139.67259 35.46523,・・・
2	神奈川県	None	None	横浜市	鶴見区	14101	POLYGON ((139.70755 35.47204, 139.70756 35.47208, 139.70719 35.47264, 139.70595 35.47454, 139.70582 35.47457, ・・・
3	神奈川県	None	None	横浜市	鶴見区	14101	POLYGON ((139.71140 35.48577, 139.71141 35.48575, 139.71108 35.48560, 139.71114 35.48551, 139.71097 35.48543, 139.71091 35.48553,・・・
4	神奈川県	None	None	横浜市	鶴見区	14101	POLYGON ((139.70710 35.47226, 139.70713 35.47228, 139.70715 35.47227, 139.70722 35.47216, 139.70722 35.47215,・・・

1275, 7

地図の表示

次に地図をプロットしてみます。 boundary.plotは境界線のみ描画します。

linewidthやedgecolorで線を調節します。

 

view plainprint?

1. gdf.boundary.plot(linewidth=0.5, edgecolor="gray", figsize=(12, 8))  
2. plt.title("神奈川県")  
3. plt.show()  

gdf.boundary.plot(linewidth=0.5, edgecolor="gray", figsize=(12, 8))
plt.title("神奈川県")
plt.show()

神奈川県の地区データが表示されました。

ValueError: numpy.dtype size changed, may indicate binary incompatibility. というエラーが出たケース

 以前動いたnumpyのコードからエラーが出るようになりました。

エラー内容

ValueError: numpy.dtype size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject

どうやら他のライブラリと互換性が無い場合に出るエラーの様で、新しくなったnumpyに対応出来ていないライブラリが使われているとのことです。

 pip freeze | grep numpy

numpy==2.0.2

ライブラリが対応するまでnumpyのバージョンを落とすなどの措置が必要です。

numpyが2.0以上なら、それ以前の最新版1.26.4に落とします。

pip uninstall numpy

pip install numpy=1.26.4

一先ず解決出来ました。

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

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

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

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

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. import geopandas as gpd  
5. import matplotlib.pyplot as plt  
6.   
7. # 世界地図のデータ  
8. world = gpd.datasets.get_path('naturalearth_lowres')  
9. # 主要都市のデータ  
10. cities = gpd.datasets.get_path('naturalearth_cities')  
11.   
12. wdf = gpd.read_file(world)  
13. cdf = gpd.read_file(cities)  
14.   
15. print(wdf.head())  
16. print()  
17. print(wdf.shape)  
18. print()  
19. print(cdf.head())  
20. print()  
21. print(cdf.shape)  

#!/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に指定することで両方プロットすることが出来ます。

view plainprint?

1. # ベースとなる世界地図  
2. bg = wdf.plot(column='continent', figsize=(18, 8))  
3. # 主要都市をマッピング  
4. cdf.plot(ax=bg, marker='o', color='firebrick')  
5.   
6. # アノテーションで主要都市の名前をマッピング  
7. for i, (city, point) in cdf.iterrows():  
8.     # pointオブジェクトから緯度経度を抽出  
9.     x1, y1, x2, y2 = point.bounds  
10.     bg.annotate(city, xy=(x1, y1))  
11.   
12. plt.show()  

# ベースとなる世界地図
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データがあるので、それらを色で段階的に分けて表示します。

人口マップ

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. import geopandas as gpd  
5. import matplotlib.pyplot as plt  
6. #日本語を表示するため  
7. import japanize_matplotlib  
8.   
9.   
10. fn = gpd.datasets.get_path('naturalearth_lowres')  
11. gdf = gpd.read_file(fn)  
12.   
13. # 人口データをint型にします  
14. gdf['population'] = gdf['pop_est'].astype("int")  
15. # 人口の多い順にソートします  
16. gdf = gdf.sort_values("pop_est", ascending=False)  
17. gdf = gdf.drop("pop_est", axis=1)  
18.   
19. print(gdf.head(10))  

#!/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

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

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

view plainprint?

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

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

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

 

各国のGDPマップ

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

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. import geopandas as gpd  
5. import matplotlib.pyplot as plt  
6. import japanize_matplotlib  
7.   
8.   
9. fn = gpd.datasets.get_path('naturalearth_lowres')  
10. gdf = gpd.read_file(fn)  
11.   
12. gdf = gdf.rename(columns={"gdp_md_est": "gdp_md"})  
13. # GDPの多い順に並べます  
14. gdf = gdf.sort_values("gdp_md", ascending=False)  
15.   
16. print(gdf.head(10))  
17.   
18. gdf.plot(column='gdp_md', cmap="jet", legend=True, figsize=(16, 8))  
19. plt.title("GDP")  
20. plt.show()  

#!/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のコロプレスマップを表示出来ました。

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

GeoPandasを使いPythonで地図を描画

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

pip install geopandas

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

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

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

 naturalearth_lowresを使用してみます。

view plainprint?

1. #!/usr/bin/env python  
2. # -*- coding: utf-8 -*-  
3.   
4. import geopandas as gpd  
5. import matplotlib.pyplot as plt  
6.   
7. # ①地図データファイルのあるパスを指定  
8. path = gpd.datasets.get_path('naturalearth_lowres')  
9.   
10. # ②ファイルからGeoDataFrameオブジェクトを生成  
11. gdf = gpd.read_file(path)  
12.   
13. print(gdf.head())  
14. print("")  
15. print(gdf.shape)  
16. print("")  
17. print(gdf.crs)  

#!/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)のことで、位置情報のルールを表します。

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

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

view plainprint?

1. gdf.plot()  
2. plt.show()  

gdf.plot()
plt.show()

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

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

view plainprint?

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

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

綺麗になりました。

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

2点間を通る直線の式

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

Pythonで書いてみます。

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

view plainprint?

1. #! /usr/bin/env python    
2.     
3. # -*- coding:utf-8 -*-    
4.      
5. #   
6. import numpy as np    
7. import matplotlib.pyplot as plt    
8. from sympy import Symbol, solve  
9.     
10. # 一次関数のデータを作る関数     
11. def linear(a, b, ran=np.arange(-10,10)):    
12.     
13.     arr = []    
14.     for x in ran:    
15.         y = x*a + b    
16.         arr.append(y)    
17.     
18.     # 式から配列データを返します  
19.     return arr    
20.     
21. ax = fig.add_subplot(1, 1, 1)    
22.      
23. # これで正方形になります   
24. ax.set_aspect('equal', adjustable='box')    
25.     
26. # 表示する幅    
27. ax.set_xlim(-10, 11)    
28. ax.set_ylim(-10, 11)    
29.     
30. # グラフの刻み    
31. ax.set_xticks(np.arange(-10, 10, 5))    
32. ax.set_yticks(np.arange(-10, 10, 5))    
33.      
34. # 補助線をonにします    
35. ax.minorticks_on()     
36.     
37. # グリッドの主線と補助線    
38. ax.grid(which="major", color="black", alpha=0.5)    
39. ax.grid(which="minor", color="gray", linestyle="--", alpha=0.5)    
40.     
41. # 数値表示を0に合わせます  
42. ax.spines['bottom'].set_position("zero")    
43. ax.spines['left'].set_position("zero")    
44.   
45. # xの生成  
46. x = np.arange(-10, 10)  
47.   
48. # 2点の位置 (-1, -5) (3, 3)  
49. x1, y1 = (-1, -5)  
50. x2, y2 = (3, 3)  
51.   
52. a = Symbol("a")  
53. b = Symbol("b")  
54.   
55. # 2点を連立方程式に変えます  
56. s1 = a * x1 + b - y1  
57. s2 = a * x2 + b - y2  
58.   
59. result = solve((s1, s2))  
60.   
61. # 得られた解  
62. print(result)  
63.   
64. # yの生成  
65. y = linear(result[a], result[b], x)  
66.   
67. # 得られた一次方程式のプロット  
68. ax.plot(x, y)  
69.   
70. # 点(-1, -5) のプロット  
71. ax.plot(x1, y1, marker="o")  
72. ax.annotate("(-1, -1)", xy=(x1, y1), color="red")  
73. # 点(3, 3) のプロット   
74. ax.plot(x2, y2, marker="o")  
75. ax.annotate("(3, 3)", xy=(x2, y2), color="red")  
76. plt.show()    

#! /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点間から一次関数のグラフが描画出来ました。