🏕️

GIS × Python Tutorial Session3 ~ shapely & pyproj~

2023/12/20に公開

Python

はじめに

この記事は「GIS × Python Tutorial」の関連記事です。
今回はPythonでGISのVectorデータを扱う際に使用する "shapely" と "pyproj" に関しての記事です。
具体的には

shapelyでのgeometry objectの作成方法
pyprojでの投影変換
shapelyの便利なメソッド

を紹介していきます。
Session5ではgeopandasを扱う予定ですが、geopandasはshapelyをベースに作成されているので、今回紹介するメソッドも使用できるはずです。

Vectorデータの種類に関しては "Session1" を参照してください。

実行環境

pythonの実行環境はGoogle Colaboratoryを想定しています。

shapelyとは

Shapelyでは、ジオメトリの操作や分析を行うことができます。二点間の距離を計測したり、エリアの面積を計測することができます。
また、wktやjsonのシリアライズ、デシリアライズができるため、地理空間データを計算するためのライブラリとして広く使われています。

pyprojとは

pyprojでは、座標系の変換や座標の計算を行うことができます。pyprojは、PROJ.4というC言語で書かれたライブラリをPythonから使えるようにしたものです。pyprojを使うことで、緯度経度座標から平面直角座標系への変換や、逆に平面直角座標系から緯度経度座標への変換ができます。

Import

import string
from typing import Iterable, List, Tuple

from matplotlib import pyplot as plt
import pyproj
import shapely
from shapely import plotting
plt.style.use('seaborn-v0_8-whitegrid')

shapely

Point Objectの作成

TupleやListで渡すも可。

geom_pt = shapely.geometry.Point(lon, lat)

MultiPoint Objectの作成

TupleやList、shapely.geometry.Pointを渡す。

geom_mpt = (
    shapely
    .geometry
    .MultiPoint([
        (lon1, lat1, alti1), 
        (lon2, lat2, alit2),
        (lon_n, lat_n, alit_n)
    ])
)

LineString Objectの作成

Listに格納した順にラインが繋がります。

geom_ln = (
    shapely
    .geometry
    .LineString([
        (lon1, lat1, alti1), 
        (lon2, lat2, alit2),
        (lon_n, lat_n, alit_n)
    ])
)

MultiLineStringの作成

geom_mln = (
    shapely
    .geometry
    .MultiLineString([
        line1, 
        line2,
        line_n
    ])
)

Polygon Objectの作成

最初の座標を最後に渡さなくても勝手に追加されます。渡しても同じです。

geom_poly = (
    shapely
    .geometry
    .Polygon([
        (lon1, lat1, alti1), 
        (lon2, lat2, alit2),
        (lon_n, lat_n, alit_n)
    ])
)

MultiPolygon Objectの作成

geom_poly = (
    shapely
    .geometry
    .MultiPolygon([
        poly1, 
        poly2
    ])
)

GeometryObjectの分解

Geometryのタイプによって操作が変わります。

MultiPointの分解

# Tuple[Tuple]での取り出し
coords = mpt.__geo_interface__.get('coordinates')

# List[shapely.geometry.Point]での取り出し
pts = list(shapely.get_parts(mpt))

LineStringの分解

# List[Tuple]での取り出し
coords = list(line.coords)

# List[shapely.geometry.Point]での取り出し
pts = [   
    shapely.get_point(line, i)
    for i in range(shapely.get_num_points(line))
]

MultiLineStringの分解

# List[List[Tuple]]での取り出し
coords = [
    list(ln.coords) 
    for ln in shapely.get_parts(mline)
]

# List[List[shapely.geometry.Point]]での取り出し
pts = []
for line in shapely.get_parts(mline):
    _pts = []
    for i in range(shapely.get_num_points(line)):
        pt = shapely.get_point(line, i)
        _pts.append(pt)
    pts.append(_pts)

Polygonの分解

# List[Tuple]での取り出し
coords = list(poly.exterior.coords)

# List[shapely.geometry.Point]での取り出し
pts = [
    shapely.get_point(poly.exterior, i)
    for i in range(shapely.get_num_points(poly.exterior))
]

MultiPolygonの分解

# List[List[Tuple]]での取り出し
coords = [
    list(poly.exterior.coords) 
    for poly in shapely.get_parts(mpoly)
]

# List[List[shapely.geometry.Point]]での取り出し
pts = []
for line in shapely.get_parts(mpoly):
    _pts = []
    for i in range(shapely.get_num_points(poly.exterior)):
        pt = shapely.get_point(poly.exterior, i)
        _pts.append(pt)
    pts.append(_pts)

pyproj

投影変換

ここでは変換器の作成と投影変換を行う関数を作成します。

投影変換する為のクラスと関数

class TransformProject(object):
    def _transformer(
        self,
        in_epsg: int,
        out_epsg: int
    ) -> pyproj.Transformer:
        # 投影変換の変換器作成
        transformer = (
            pyproj
            .Transformer
            .from_crs(
                f'epsg:{in_epsg}', 
                f'epsg:{out_epsg}', 
                always_xy=True
            )
        )
        return transformer

    def _x_y_to_xy(
        self,
        lons: Iterable[float],
        lats: Iterable[float]
    ) -> List[Tuple[float]]:
        # 経度のListと緯度のListから経緯度が入ったTupleのListを作成する
        return [
            (lon, lat) 
            for lon, lat in zip(lons, lats)
        ]

    def _xy_to_x_y(
        self,
        coords_lst: Iterable[float] | Iterable[Iterable[float]]
    ) -> List[List[float]]:
        # 経緯度が入ったTupleのListから経度のListと緯度のListを作成する
        if isinstance(coords_lst[0], Iterable):
            lons = [coords[0] for coords in coords_lst]
            lats = [coords[1] for coords in coords_lst]
        else:
            lons = coords_lst[0]
            lats = coords_lst[1]
        return [lons, lats]


def transform_coords(
    in_epsg: int,
    out_epsg: int,
    lon: float | Iterable[float]=None,
    lat: float | Iterable[float]=None,
    coords: Iterable[Iterable[float]]=None,
    return_type_x_y: bool=True
) -> Iterable[Iterable[float]]:
    """
    in_epsg(int): 入力座標のEPSGコード
    out_epsg(int): 出力座標のEPSGコード
    lon(float | Iterable[float]): 経度(x)、経度(x)のList
    lat(float | Iterable[float]): 緯度(y)、緯度(y)のList
    coords(Iterable[Iterable[float]]): 経緯度のTupleが入ったList
    return_type_x_y(bool): Trueならば経度と緯度を別なListで返す
                           Falseならば経緯度のTupleが入ったListを返す
    """
    tp = TransformProject()
    transformer = tp._transformer(in_epsg, out_epsg)
    if coords:
        lon, lat = tp._xy_to_x_y(coords)
    # 投影変換
    coords = transformer.transform(lon, lat)
    if return_type_x_y:
        # List[lon], List[lat]で返す
        return coords
    else:
        if isinstance(coords[0], Iterable):
            # List[tuple(lon, lat)]で返す
            return tp._x_y_to_xy(coords[0], coords[1])
        else:
            return coords

lon = 140.7418870
lat = 40.82524732
lons = [lon, lon+0.01]
lats = [lat, lat+0.01]
projected_coords = transform_coords(4326, 6678, lons, lats, return_type_x_y=False)
print(projected_coords)

output
[
(-7713.368225111149, 91632.44828291438),
(-6868.850338328306, 92742.01169168559)
]

shapelyの便利なメソッド

ここからはshapelyの便利なメソッドを紹介していきます。

データの準備

まずは距離の計算を簡単にする為に平面直角座標系に投影変換します。

x, y = transform_coords(4326, 6678, lon, lat)
point_a = shapely.geometry.Point(x, y)
point_b = shapely.geometry.Point(x + 100, y)
point_c = shapely.geometry.Point(x + 50, y - 50)
point_d = shapely.geometry.Point(x + 50, y + 50)
points = [point_a, point_b, point_c, point_d]

# ここはあまり関係ない、単なるラベル作成用の関数
def alhpabet_lst(num) -> str:
    # アルファベットの生成
    alphabet = list(string.ascii_uppercase)
    return alphabet[num]

Geometryのplot

shapelyには "shapely.plotting" というplot用のメソッドがあります。
matplotlibを使用した事がある方であればドキュメントを読みに行かなくても使用できるでしょう。

# Pointの位置をPlot
for i, pt in enumerate(points):
    plotting.plot_points(pt, label=f'point-{alhpabet_lst(i)}')

plt.legend()
plt.show()
# ラベルなしでもいいのであれば
# plotting.plot_points(points)でも可

バッファーの作成 [geom.buffer]

bufferは各ジオメトリオブジェクトで作成可能です。

geom.buffer(
	# Bufferの距離を指定（平面直角座標系なのでメートルで指定）
	distance=10,
	# 解像度を指定（数値を大きくすれば解像度が上がる）
	resolustion=1,
	# bufferの形状を指定[round | square | flat]
	cap_style='round',
	# Bufferを片側だけに作成する（distance引数が正の距離は左、負の距離は右に）
	# ドキュメントだとsingle_sideになっている
	single_sided=False,
)

適当にBufferを作成

buffers = []
for pt in points:
    buff = pt.buffer(60)
    buffers.append(buff)

multi_buffer = shapely.geometry.MultiPolygon(buffers)
# Plot
plotting.plot_polygon(
    multi_buffer, 
    add_points=False, 
    color='gray', 
    label='Buffer'
)
plotting.plot_points(
    geom=points, 
    label='Points'
)
plt.legend()
plt.show()

resolutionを試す

cap_styleを試す

Lineの端からBufferを伸ばしたくない場合は"flat"が使えそうですね。

single_sidedを試す

single_sidedをTrueにすると、cap_styleが強制的に"flat"になります。

LineStringを作成する際にPointsを繋ぐ順番を入れ替えると、左右の意味が変わるので注意して下さい。

line = (
    shapely
    .geometry
    .LineString(
        points[: 3]
    )
)

line_reverse = (
    shapely
    .geometry
    .LineString(
        list(reversed(points[: 3]))
    )
)
# 両方とも左側にBufferを作成する
l_buff = line.buffer(10, single_sided=True)
lr_buff = line_reverse.buffer(10, single_sided=True)

fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
plotting.plot_line(line, ax[0], add_points=False)
plotting.plot_line(line_reverse, ax[1], add_points=False)
plotting.plot_polygon(l_buff, ax[0], add_points=False)
plotting.plot_polygon(lr_buff, ax[1], add_points=False)

Bufferは負の値を指定する事も出来ます。

fig = plt.figure(figsize=(10, 10))
plotting.plot_polygon(
    buffers[0], 
    add_points=False, 
    fill=None,
    label='Original'
)

# 内側に小さなBufferを作成
in_buff = buffers[0].buffer(-10)
plotting.plot_polygon(
    in_buff, 
    add_points=False, 
    color='red',
    label='Small'
)
plt.legend()
plt.show()

差分の抽出 [geom.difference]

geometry同士の差分を計算し、差分のgeometryを返します。
このメソッドを使用して内向きのBufferを作成してみましょう。

geom1.difference(geom2)

out_buff = buffers[0]
in_buff = out_buff.buffer(-10)
inner_buff = out_buff.difference(in_buff)

fig = plt.figure(figsize=(8, 8))
plotting.plot_polygon(
    inner_buff,
    add_points=False,
    color='red',
    label='Innter buffer'
)
plotting.plot_polygon(
    out_buff, 
    add_points=False, 
    fill=None,
    lw=3,
    label='Original'
)
plt.legend()
plt.show()

重なりの抽出 [geom.intersection]

intersectionでは重なり合う面のgeometryを返します。
作成済みの4つのbufferの重なりを求めてみましょう。

geom1.intersection(geom2)

# 重なり合う面を抽出する
intersect = buffers[0]
for buff in buffers[1: ]:
    intersect = intersect.intersection(buff)


# Plot
fig = plt.figure(figsize=(8, 8))
plotting.plot_polygon(
    shapely.geometry.MultiPolygon(buffers),
    add_points=False,
    color='gray',
    label='buffers'
)
plotting.plot_polygon(
    intersect,
    add_points=False,
    color='red',
    label='buffer intersection'
)
plotting.plot_points(
    geom=points,
    label='Points'
)
plt.legend()
plt.show()

長さの計測 [geom.length]

line.length

intersect.length

面積の計算 [geom.area]

geom.area

xyの最小値と最大値を計算 [geom.bounds]

geom.bounds

Geometryのxyの最小値と最大値からBoxPolygonを作成 [shapely.box]

box_poly = shapely.box(*intersect.bounds)
plotting.plot_polygon(box_poly, color='red', label='box')
plotting.plot_polygon(intersect, add_points=False, fill=None, label='intersect')
plt.legend();

中心を計算 [geom.centroid]

geom.centroid

Geometry同士のユークリッド距離計算 [geom1.distance(geom2) | shapely.distance(geom1, geom2)]

shapelyではユークリッド距離の他にハウスドルフ距離やフレチェット距離も計算できますが、解説すると長くなるので別なSessionで解説する予定です。

# データの準備
points_a = shapely.geometry.MultiPoint([
    shapely.geometry.Point(0, 50),
    shapely.geometry.Point(-10, 30)
])
points_b = shapely.geometry.MultiPoint([
    shapely.geometry.Point(20, 50),
    shapely.geometry.Point(40, 80),
    shapely.geometry.Point(60, 0),
    shapely.geometry.Point(-5, 20)
])

# Plot
for i, p in enumerate(shapely.get_parts(points_a)):
    plotting.plot_points(p, label='A', color='red')
    plt.text(p.x + 3, p.y, f"$A_{i}$")

for i, p in enumerate(shapely.get_parts(points_b)):
    plotting.plot_points(p, label='B', color='blue')
    plt.text(p.x + 3, p.y, f"$B_{i}$")

plt.axis('equal')
plt.show()

for i, a in enumerate(shapely.get_parts(points_a)):
    for j, b in enumerate(shapely.get_parts(points_b)):
        euclidean = a.distance(b)
        print(f"A{i} から B{j}までのユークリッド距離は: {round(euclidean, 2)}")

output
A0 から B0までのユークリッド距離は: 20.0
A0 から B1までのユークリッド距離は: 50.0
A0 から B2までのユークリッド距離は: 78.1
A0 から B3までのユークリッド距離は: 30.41
A1 から B0までのユークリッド距離は: 36.06
A1 から B1までのユークリッド距離は: 70.71
A1 から B2までのユークリッド距離は: 76.16
A1 から B3までのユークリッド距離は: 11.18

Geometryの結合 [geom1.union(geom2)]

result = inner_buff.union(buffers[1])
plotting.plot_polygon(result, add_points=False)

geojson formatのDictを作成 [geom.geo_interface]

geom = point_a.buffer(10, resolution=1)
geom.__geo_interface__

{
    'type': 'Polygon',
    'coordinates': (
        (
            (-7703.368225111149, 91632.44828291438),
            (-7713.368225111149, 91622.44828291438),
            (-7723.368225111149, 91632.44828291438),
            (-7713.368225111149, 91642.44828291438),
            (-7703.368225111149, 91632.44828291438)
        ),
    )
}

Geometryの回転 [shapely.affinity.rotate]

shapely.affinity.rotate(geom, d)

# オブジェクトの回転
fig = plt.figure(figsize=(7, 7))
plotting.plot_polygon(result, add_points=False, label='original')
rotate = shapely.affinity.rotate(result, 90)
plotting.plot_polygon(rotate, add_points=False, color='red', label='rotate')
plt.legend()
plt.show()

Geometryのスケーリング [shapely.affinity.scale]

# オブジェクトの縮小
fig = plt.figure(figsize=(7, 7))
plotting.plot_polygon(result, add_points=False, label='original')
reduction = shapely.affinity.scale(result, xfact=0.75, yfact=0.75)
plotting.plot_polygon(reduction, add_points=False, color='red', label='reduction')
plt.legend()
plt.axis('equal')
plt.show()

# オブジェクトの拡大
fig = plt.figure(figsize=(7, 7))
plotting.plot_polygon(result, add_points=False, label='original')
expansion = shapely.affinity.scale(result, xfact=1.1, yfact=1.1)
plotting.plot_polygon(expansion, add_points=False, color='red', label='expansion')
plt.legend()
plt.axis('equal')
plt.show()

Geometryの移動 [shapely.affinity.translate]

# オブジェクトの移動
fig = plt.figure(figsize=(7, 7))
plotting.plot_polygon(result, add_points=False, label='original')
shift = shapely.affinity.translate(result, xoff=20, yoff=10)
plotting.plot_polygon(shift, add_points=False, color='red', label='shift')
plt.legend()
plt.axis('equal')
plt.show()

おわりに

いかがでしたでしょうか。shapelyのメソッドはまだまだ沢山ありますが、shapely単体で使用するのはこれくらいかな?というのを紹介しました。
今後のSessionで登場するgeopandasのgeometryではshapelyが使用されているので、これらのメソッドを使用する事も可能です。
次回のSessionではgeopandasについての記事を執筆する予定です。ちょっと続けて投稿したので期間が開くかもしれません。

はじめに

実行環境

shapelyとは

pyprojとは

Import

shapely

Point Objectの作成

MultiPoint Objectの作成

LineString Objectの作成

MultiLineStringの作成

Polygon Objectの作成

MultiPolygon Objectの作成

GeometryObjectの分解

MultiPointの分解

LineStringの分解

MultiLineStringの分解

Polygonの分解

MultiPolygonの分解

pyproj

投影変換

shapelyの便利なメソッド

データの準備

Geometryのplot

バッファーの作成 [geom.buffer]

適当にBufferを作成

resolutionを試す

cap_styleを試す

single_sidedを試す

Bufferは負の値を指定する事も出来ます。

差分の抽出 [geom.difference]

重なりの抽出 [geom.intersection]

長さの計測 [geom.length]

面積の計算 [geom.area]

xyの最小値と最大値を計算 [geom.bounds]

Geometryのxyの最小値と最大値からBoxPolygonを作成 [shapely.box]

中心を計算 [geom.centroid]

Geometry同士のユークリッド距離計算 [geom1.distance(geom2) | shapely.distance(geom1, geom2)]

Geometryの結合 [geom1.union(geom2)]

geojson formatのDictを作成 [geom.geo_interface]

Geometryの回転 [shapely.affinity.rotate]

Geometryのスケーリング [shapely.affinity.scale]

Geometryの移動 [shapely.affinity.translate]

おわりに

Discussion