🏕️
GIS × Python Tutorial 5.1 ~ geopandas ことはじめ GeoDataFrameの基本 ~
はじめに
この記事は「GIS × Python Tutorial」の関連記事です。
今回はgeopandas.GeoDataFrameの基本について解説します。pythonを使用する方であればpandasを使用した事がある方が多いかと思いますので、分からない部分はgeopandasドキュメントを見て行けば躓かずに理解できるかと思います。またgeopandasのgeometryはshapelyのgeometryオブジェクトが入力されるので、個別のメソッドはshapelyの公式ドキュメント、あるいは前回の記事を参考にしてください。
geopandas とは
geopandasは、pandasの拡張であり地理空間情報を扱いやすくした DataFrame を操作する事が可能です。ドキュメントも非常に読みやすいと思いますので、このシリーズを読んでもっと詳しく知りたい方はドキュメントを読んでみてください。またチュートリアルも充実しているようなので GitHub のリンクも貼っておきます。また今回の実行環境やこの記事を執筆する際に使用したNotebookはGitHubに公開しているので、動かしてみたい方はどうぞ
Install geopandas
コード実行の準備
Import
cell 1
# 今回使用するライブラリのインポート
import fiona
import geopandas as gpd
from matplotlib import pyplot as plt
import numpy as np
import shapely
plt.style.use('ggplot')
データの作成
適当にデータを作成します。
cell 2
# 弘前城の位置
x1 = 140.46362773837438
y1 = 40.60790364299233
point1 = shapely.Point(x1, y1)
# 弘前城二の丸の位置
x2 = 140.46539395937307
y2 = 40.60780032475679
point2 = shapely.Point(x2, y2)
# ランダムな点群を作成
length = 200
alphabet = 'ABCDE'
xs = np.random.normal(x1, abs(x1 - x2), length)
ys = np.random.normal(y1, abs(y1 - y2), length)
points = [shapely.Point(_x, _y) for _x, _y in zip(xs, ys)]
# アルファベットからランダムに値を取得
indexes = np.random.randint(0, len(alphabet), length)
codes = [alphabet[idx] for idx in indexes]
GeoDataFrameの作成
geopandas.GeoDataFrameの作成はほとんどpandas.DataFrameと同じです。違うのは geometry と crs を設定する点です。geometry には shapely.geometry.XXX を List などに入れて渡します。crs には EPSG コードや WKT-CRS文字列等の座標系が識別できるものを渡します。
cell 3
IN_EPSG = 'EPSG:4326'
gdf = (
gpd.GeoDataFrame(
data={'CODE': codes},
# 今回はgeometryにList[shapely.Point]を渡す。
geometry=points,
crs=IN_EPSG
)
)
print(f"shape: {gdf.shape}")
print(gdf.head().to_markdown())
shape: (200, 2)
| CODE | geometry | |
|---|---|---|
| 0 | E | POINT (140.45955045678622 40.60782384942577) |
| 1 | C | POINT (140.4629219485244 40.60792626484476) |
| 2 | C | POINT (140.4637344764471 40.60783463065041) |
| 3 | A | POINT (140.46080018253036 40.607916728840046) |
| 4 | E | POINT (140.45998875601498 40.60794931474045) |
geometry
geometry の列を見てみましょう。geometry の列はgeopandas.geoseries.GeoSeriesになっているのが確認できます。この状態の時には geopadnas は地理的空間的なメソッドを扱う事が出来る様になっています。pandas のメソッドを扱ったりするとこれが外れたりするので、エラーが出たら確認してみましょう。
cell 4
print(f"""
<< dtypes >>
{gdf.dtypes}
<< series type >>
CODE = {type(gdf['CODE'])}
geometry = {type(gdf.geometry)}
""")
output 4
<< dtypes >>
CODE object
geometry geometry
dtype: object
<< series type >>
CODE = <class 'pandas.core.series.Series'>
geometry = <class 'geopandas.geoseries.GeoSeries'>
geometry が設定されていない場合は以下の様にエラーが出ます。
cell 5
_gdf = gpd.GeoDataFrame(data={'CODE': codes, 'SHAPE': points})
_gdf.geometry
output 5
File X:\XX\envs\lidar\lib\site-packages\geopandas\geodataframe.py:236, in GeoDataFrame._get_geometry(self)
229 else:
230 msg += (
231 "\nThere are no existing columns with geometry data type. You can "
232 "add a geometry column as the active geometry column with "
233 "df.set_geometry. "
234 )
--> 236 raise AttributeError(msg)
237 return self[self._geometry_column_name]
AttributeError: You are calling a geospatial method on the GeoDataFrame, but the active geometry column to use has not been set.
There are no existing columns with geometry data type. You can add a geometry column as the active geometry column with df.set_geometry.
geometry を設定してみましょう。実は geometry の列名は何でも構いません。
cell 6
_gdf.set_geometry('SHAPE', inplace=True)
print(f"""
<< dtypes >>
{_gdf.dtypes}
<< series type >>
CODE = {type(_gdf['CODE'])}
geometry = {type(_gdf.geometry)}
""")
output 6
<< dtypes >>
CODE object
SHAPE geometry
dtype: object
<< series type >>
CODE = <class 'pandas.core.series.Series'>
geometry = <class 'geopandas.geoseries.GeoSeries'>
CRSの確認
crs を設定していればGeodataFrame.crsのメソッドでCRSを確認する事が出来ます。
cell 7
gdf.crs
output 7
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
wkt-crs も見てみます。
cell 8
print(gdf.crs.to_wkt(pretty=True))
output 8
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
EPSGコードも出力してみましょう。
cell 9
print(gdf.crs.to_epsg())
output 9
4326
CRSが設定されていない場合
CRSが設定されていない場合は設定する事が可能です。
cell 10
gdf.crs = None
print(f"EPSG 1: {gdf.crs}")
gdf.set_crs(IN_EPSG, inplace=True)
print(f"EPSG 2: {gdf.crs.to_epsg()}")
output 10
EPSG 1: None
EPSG 2: 4326
投影変換
geopandas での投影変換は非常に簡単です。しかし、CRS が設定されていない場合はエラーになるので、上で書いた様にCRSは確認する様にしましょう。
cell 11
OUT_EPSG = 'EPSG:6678'
gdf_jgd = gdf.to_crs(OUT_EPSG)
print(f"""
IN_EPSG: {gdf.crs.to_epsg()}
OUT_EPSG: {gdf_jgd.crs.to_epsg()}
""")
output 11
IN_EPSG: 4326
OUT_EPSG: 6678
UTM座標系の推定
estimate_utm_crsメソッドで UTM 座標系を推定する事が出来ます。都道府県を跨ぐような中規模のデータセットを使用する場合には非常に便利です。
estimate_utm_crsは "datum_name" を指定する事が出来ます。デフォルトは 'WGS 84' ですがいくつか見てみましょう。
cell 12
gdf.estimate_utm_crs()
output 12
<Projected CRS: EPSG:32654>
Name: WGS 84 / UTM zone 54N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 138°E and 144°E, northern hemisphere between equator and 84°N, onshore and offshore. Japan. Russian Federation.
- bounds: (138.0, 0.0, 144.0, 84.0)
Coordinate Operation:
- name: UTM zone 54N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
cell 13
gdf.estimate_utm_crs(datum_name='JGD2000')
output 13
<Projected CRS: EPSG:3100>
Name: JGD2000 / UTM zone 54N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Japan - between 138°E and 144°E, onshore and offshore.
- bounds: (138.0, 17.63, 144.0, 46.05)
Coordinate Operation:
- name: UTM zone 54N
- method: Transverse Mercator
Datum: Japanese Geodetic Datum 2000
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
cell 14
gdf.estimate_utm_crs(datum_name='JGD2011')
output 14
<Projected CRS: EPSG:6691>
Name: JGD2011 / UTM zone 54N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Japan - between 138°E and 144°E, onshore and offshore.
- bounds: (138.0, 17.63, 144.0, 46.05)
Coordinate Operation:
- name: UTM zone 54N
- method: Transverse Mercator
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
これを使用して投影変換してみましょう。
cell 15
wkt_crs = gdf.estimate_utm_crs(datum_name='JGD2011')
gdf_utm = gdf.to_crs(wkt_crs)
print(f"EPSG: {gdf_utm.crs.to_epsg()}")
output 15
EPSG: 6691
shapely のメソッド
GeoSeries では shapely のメソッドを使用する事が出来ます。これが geopandas が強力なライブラリーである 1つの理由です。
cell 16
# 各Pointから10mのBufferを作成
buffers = gdf_jgd.geometry.buffer(10)
cell 17
# 弘前城の位置を入力したGeoDataFrameを作成
base_gdf_jgd = (
gpd.GeoDataFrame(
data={'CODE': ['BASE']},
geometry=[point1],
crs=IN_EPSG
)
.to_crs(OUT_EPSG)
)
# 各Pointと弘前城との距離を計測
gdf_jgd['from_base'] = gdf_jgd.distance(base_gdf_jgd.geometry.iloc[0])
print(gdf_jgd.head().to_markdown())
| CODE | geometry | from_base | |
|---|---|---|---|
| 0 | E | POINT (-31630.825010381694 67553.45817584553) | 345.151 |
| 1 | C | POINT (-31345.46902112123 67563.62391520104) | 59.7798 |
| 2 | C | POINT (-31276.752721893547 67553.16024279944) | 11.8452 |
| 3 | A | POINT (-31525.02493539584 67563.3227825572) | 239.284 |
| 4 | E | POINT (-31593.675417477116 67567.23188986423) | 307.988 |
データの出力
今回は geojson で出力する場合と geopackege で出力する方法の 2種類を見てみましょう。
-
geojson はファイル名のみ指定すれば保存できます。
-
geopackage はファイル名の他に Layer 名を指定する必要があります。
cell 18
# geojsonで出力する場合
OUT_FILE_GEOJSON = r'../datasets/session5/Hirosaki.geojson'
gdf.to_file(OUT_FILE_GEOJSON, driver='GeoJSON')
# geopackegeに出力する場合
OUT_FILE_GEOPACKEGE = r'../datasets/session5/Hirosaki.gpkg'
LAYER1 = 'RANDOM_DATA'
gdf.to_file(OUT_FILE_GEOPACKEGE, layer=LAYER1, driver='GPKG')
LAYER2 = 'BasePoints'
gdf.to_file(OUT_FILE_GEOPACKEGE, layer=LAYER2, driver='GPKG')
データの読み込み
geojson の読み込みは非常に簡単ですが、geopackage は Layer 名を指定する必要があります。
cell 19
_ = gpd.read_file(OUT_FILE_GEOJSON)
print(_.head().to_markdown())
| CODE | geometry | |
|---|---|---|
| 0 | E | POINT (140.45955045678622 40.60782384942577) |
| 1 | C | POINT (140.4629219485244 40.60792626484476) |
| 2 | C | POINT (140.4637344764471 40.60783463065041) |
| 3 | A | POINT (140.46080018253036 40.607916728840046) |
| 4 | E | POINT (140.45998875601498 40.60794931474045) |
cell 20
# fionaでgeopackageに保存されているLayer名を取得
layers = fiona.listlayers(OUT_FILE_GEOPACKEGE)
print(f"Layers: {layers}")
cell 20
Layers: ['RANDOM_DATA', 'BasePoints']
cell 21
gdf = gpd.read_file(OUT_FILE_GEOPACKEGE, layer=layers[0])
print(gdf.head().to_markdown())
| CODE | geometry | |
|---|---|---|
| 0 | E | POINT (140.45955045678622 40.60782384942577) |
| 1 | C | POINT (140.4629219485244 40.60792626484476) |
| 2 | C | POINT (140.4637344764471 40.60783463065041) |
| 3 | A | POINT (140.46080018253036 40.607916728840046) |
| 4 | E | POINT (140.45998875601498 40.60794931474045) |
おわりに
今回はGeoDataFrameの作成と投影法の定義など基本の操作について紹介しました。次回のチュートリアル 5.2 では空間検索について書きたいと思います。
参考
Discussion