はじめに

この記事は「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

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