PythonでのGIS・位置情報データのハンドリング実践
ShapelyやGeoPandas、Folium、Plotlyなどによるさまざまなデータ処理・可視化手法を学ぶ
OSMNXはopenstreetmap による最短経路の算出
内容の実践はdocker を使ってるけど
Pipでインストールも可能
Esri japanでGISの活用事例が載っている
今回はpythonでやるのでQGISは使わない
シェープファイルのサイズはそれぞれ２GBまで
フィールド名は英数字１０文字
日本語なら５文字まで
geopandas関連
pip install geopandas
でインストール
次に地図関連
pip install folium pip install matplotlib
で可視化ツールのインストール
pip install jupyter
でインストール
jupyter-notebook
で起動
以下コード
import geopandas as gpd import folium import matplotlib.pyplot as plt
でインポート
path = gpd.datasets.get_path('nybb')
gdf = gpd.read_file(path)
gdf.head()
でデータ取得して表示するはずが
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[3], line 1
----> 1 path = gpd.datasets.get_path('nybb')
2 gdf = gpd.read_file(path)
3 gdf.head()
File ~/.pyenv/versions/3.11.0/lib/python3.11/site-packages/geopandas/datasets/__init__.py:18, in get_path(dataset)
12 error_msg = (
13 "The geopandas.dataset has been deprecated and was removed in GeoPandas "
14 f"1.0. You can get the original '{dataset}' data from "
15 f"{ne_message if 'natural' in dataset else nybb_message}"
16 )
17 if dataset in _prev_available:
---> 18 raise AttributeError(error_msg)
19 else:
20 error_msg = (
21 "The geopandas.dataset has been deprecated and "
22 "was removed in GeoPandas 1.0. New sample datasets are now available "
23 "in the geodatasets package (https://geodatasets.readthedocs.io/en/latest/)"
24 )
AttributeError: The geopandas.dataset has been deprecated and was removed in GeoPandas 1.0. You can get the original 'nybb' data from the geodatasets package.
from geodatasets import get_path
path_to_file = get_path('nybb')
原因は
GeoPandas 1.0 以降では gpd.datasets.get_path(‘nybb’) が削除された ため
代わりに、geodatasets パッケージを使う必要がある
pip install geodatasets
でインストール
インポート文を
import geopandas as gpd import folium import matplotlib.pyplot as plt from geodatasets import get_path
にして
# path = gpd.datasets.get_path('nybb')
# データセットの取得
path = get_path("nybb")
gdf = gpd.read_file(path)
gdf.head()
これで
Downloading file 'nybb_16a.zip' from 'https://www.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_16a.zip' to '/Users/snowpool/Library/Caches/geodatasets'. Extracting 'nybb_16a/nybb.shp' from '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip' to '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip.unzip' Extracting 'nybb_16a/nybb.shx' from '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip' to '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip.unzip' Extracting 'nybb_16a/nybb.dbf' from '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip' to '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip.unzip' Extracting 'nybb_16a/nybb.prj' from '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip' to '/Users/snowpool/Library/Caches/geodatasets/nybb_16a.zip.unzip' BoroCode BoroName Shape_Leng Shape_Area geometry 0 5 Staten Island 330470.010332 1.623820e+09 MULTIPOLYGON (((970217.022 145643.332, 970227.... 1 4 Queens 896344.047763 3.045213e+09 MULTIPOLYGON (((1029606.077 156073.814, 102957... 2 3 Brooklyn 741080.523166 1.937479e+09 MULTIPOLYGON (((1021176.479 151374.797, 102100... 3 1 Manhattan 359299.096471 6.364715e+08 MULTIPOLYGON (((981219.056 188655.316, 980940.... 4 2 Bronx 464392.991824 1.186925e+09 MULTIPOLYGON (((1012821.806 229228.265, 101278...
というように無事に表示される
次に表ではなく図で描画
gdf.plot()
で簡単に描画できる
次に座標情報を表示
gdf.crs
次に座標系の変換
gdf = gdf.to_crs(epsg=4326) print(gdf.crs) gdf.head()
これで座標が緯度経度に変わっている
EPSG:4326 BoroCode BoroName Shape_Leng Shape_Area geometry 0 5 Staten Island 330470.010332 1.623820e+09 MULTIPOLYGON (((-74.05051 40.56642, -74.05047 ... 1 4 Queens 896344.047763 3.045213e+09 MULTIPOLYGON (((-73.83668 40.59495, -73.83678 ... 2 3 Brooklyn 741080.523166 1.937479e+09 MULTIPOLYGON (((-73.86706 40.58209, -73.86769 ... 3 1 Manhattan 359299.096471 6.364715e+08 MULTIPOLYGON (((-74.01093 40.68449, -74.01193 ... 4 2 Bronx 464392.991824 1.186925e+09 MULTIPOLYGON (((-73.89681 40.79581, -73.89694 ...
緯度経度に変換すれば foliumで地図表示ができる
まずは地図の表示
m = folium.Map(location=[40.7,-73.96],zoom_start=10,tiles='CaltoDB positiuon') m
としたらエラー
単純に
CartoDB positron
の間違い
m = folium.Map(location=[40.7,-73.96],zoom_start=10,tiles='CartoDB positron') m
が正解
次に地図上にgeopandas情報を描画する
まずはこの情報を上から見ていくので
for __, r in gdf.iterrows():
sim_geo = gpd.GeoSeries(r['geometory']).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['BoroName']).add_to(geo_j)
geo_j.add_to(m)
m
今回もエラー
原因は今回もスペルミス
geometory
じゃなくて
geometry
が正解
これで地図の上に情報が表示される
geopandas関連その２
今回は shapely を使う
import pandas as pd import geopandas as gpd import matplotlib as plt from shapely.geometry import Point,LineString, Polygon
まずはPoint の練習
#Point point = Point(0.0,1.0) point
これで座標を指定するとそこに⭕️が表示される
この座標を出すなら
print(point)
とすれば
POINT (0 1)
と座標を表示できる
ジオメトリタイプを調べるなら
point.geom_type
というように
.geom_type
をつけることでタイプの判別ができる
次にLineStringの練習
line = LineString([(0,0),(1,1),(3,0),(5,2)]) line
これで折れ線が表示される
print(line)
で
LINESTRING (0 0, 1 1, 3 0, 5 2)
というように配列の座標を見れる
次はPolygon
polygon = Polygon([(0,0),(1,1),(3,0)]) polygon
これで三角形が表示される
polygon = Polygon([(0,0),(0,1),(1,1),(1,0)]) polygon
これで四角形
これらをまとめて行うのがジオメトリコレクション
#multipoint from shapely.geometry import MultiPoint,MultiLineString,MultiPolygon
で必要なインポートを行う
points = MultiPoint([(0,0),(1,0)]) points
で２つの座標に⭕️が表示される
coords = [((0,0),(1,1)),((-1,0),(1,0))] lines = MultiLineString(coords) lines
これで２つの線が憑依される
a = Polygon([(0,0),(0,1),(1,1),(1,0)]) b =Polygon([(2,0),(2,1),(3,1),(3,0)]) polygons = MultiPolygon([a,b]) polygons
で２つの四角が表示される
次に LineString の長さ
l = LineString([(1,-1),(1,0),(2,0),(2,1)]) l
この長さを求めるには
.length を使う
これは道路などで使うことになる
次に polygonの面積
p=Polygon([(3,-1),(4,0),(3,1)]) p
で三角形
面積は
areaで求めることができるので
p.area
とすれば
1.0
となる
次にbuffer
これは Point LineString に幅を与えてpolygonにするもの
s = gpd.GeoSeries(
[
Point(0,0),
LineString([(1,-1),(1,0),(2,0),(2,1)]),
Polygon([(3,-1),(4,0),(3,1)])
])
s
結果は
0 POINT (0 0) 1 LINESTRING (1 -1, 1 0, 2 0, 2 1) 2 POLYGON ((3 -1, 4 0, 3 1, 3 -1)) dtype: geometry
次に
s.buffer(0.5)
とすると
0 POLYGON ((0.5 0, 0.49759 -0.04901, 0.49039 -0.... 1 POLYGON ((0.5 0, 0.50241 0.04901, 0.50961 0.09... 2 POLYGON ((2.5 -1, 2.5 1, 2.50241 1.04901, 2.50... dtype: geometry s.buffer(0.5)[0]
で円ができる
これは半径0.5の円
次にLineString での buffer
s.buffer(0.5)[1]
これで線の周りが追加されて図みたいになっている
この時に角はまるくなる
同様にpolygonも試す
s.buffer(0.5)[2]
これで三角形の角がまるくなる
これらは後々、ジオメトリの空間結合で使う
例えば予測地点の１００m以内のみ表示などで使う