geopandasでshpファイルのポリゴンを構成しているポイント座標（ジオメトリ）を抽出したい

Question

geopandasでshpファイルのポリゴンを構成しているポイント座標（ジオメトリ）を抽出したい

農地データと人工衛星のマイクロ波の照射方向の関係性について研究しています。
geopandasによってshpファイル（ポリゴン）を抽出しました。
そこでジオメトリのポイント座標を演算のために抽出しようと考えています。

###実行環境
実行環境はWindows10、Annaconda（Spyder）を使用しています。

該当のソースコード

import geopandas as gpd

#shpファイルの読み込み
yoshi_2019=gpd.read_file("yoshikawa_2019.shp")

yoshi_points = yoshi_2019['geometry']
print(yoshi_points)

出力結果

上記のコードを実行すると、以下の結果が得られました。

0       POLYGON Z ((273474.470 4123032.881 0.000, 2734...
1       POLYGON Z ((266669.419 4122033.251 0.000, 2666...
2       POLYGON Z ((266698.965 4122029.388 0.000, 2667...
3       POLYGON Z ((266659.482 4121927.719 0.000, 2666...
4       POLYGON Z ((266626.999 4121931.518 0.000, 2666...
                       
3693    POLYGON Z ((266500.679 4121835.558 0.000, 2665...
3694    POLYGON Z ((266524.600 4121831.832 0.000, 2665...
3695    POLYGON Z ((266686.077 4121924.596 0.000, 2667...
3696                                                 None
3697                                                 None

Name: geometry, Length: 3698, dtype: geometry

実現させたいこと

私はPOLIGON　Z　（）のポリゴンを構成している座標を抽出したいです。
おそらく　（273474.470 4123032.881 0.000）を取り出すことが出来れば
、演算できるのではないかと考えています。
もし分かる方がおられましたら、教えて頂きたいです。
よろしくお願いします。

Answer 1

"POLYGON Z ((273474.470 4123032.881 0.000, 2734..." と表示されているデータの実体は、GeoPandasが依存するShapelyライブラリのオブジェクトです。
ここでは格納されているのがポリゴンデータですので、Polygonオブジェクトです。

>>> yoshi_2019['geometry'][0]
<shapely.geometry.polygon.Polygon object at 0x7fc870d37950>
>>> yoshi_2019['geometry'][0].geom_type
Polygon

Shapelyはジオメトリの操作と分析のためのライブラリです。Shapelyライブラリのオブジェクトに用意されているメソッドや属性を使って、ジオメトリのポイント座標を抽出することができます。

Shapelyのジオメトリオブジェクトにはinterior、boundary、exteriorの3つの属性があり、さらにそれぞれのcoords属性によって座標値を取得できます。
単純なポリゴンの場合は、exteriorから取得すれば良いでしょう。

# サンプルデータの作成
>>> from shapely.geometry import Polygon
>>> p = Polygon([(0, 0, 1), (1, 1, 2), (1, 0, 2), (0, 0, 1)])
>>> p
<shapely.geometry.polygon.Polygon object at 0x7fc87064d1d0>
>>> print(p)
POLYGON Z ((0 0 1, 1 1 2, 1 0 2, 0 0 1))

# exteriorからポイント座標をリストとして抽出
>>> list(p.exterior.coords)
[(0.0, 0.0, 1.0), (1.0, 1.0, 2.0), (1.0, 0.0, 2.0), (0.0, 0.0, 1.0)]

どのような形で座標を抽出するのが良いのかは、質問者さんのやりたいことに依存します。たとえばX座標、Y座標、Z座標ごとにまとめて取得したいときはlist(zip(*p.exterior.coords)) とすることで可能です。

単純に元のGeoPandasデータフレームに追加する場合は次のようにすることができます。'coordinates’という新しい列を作ってそこに追加してみました。

>>> yoshi_2019['coordinates'] = yoshi_2019['geometry'].apply(lambda p:list(p.exterior.coords))

ポリゴンデータが内部に穴を持っている場合は、内部の穴のポリゴンの座標値も取得するのかとか（この値はinteriorsに格納されています）、マルチポリゴンだった場合はそれぞれのポリゴンごとに値を取得しないといけないなど、それに応じた考慮が必要になります。

Shapely — Shapely documentation
https://shapely.readthedocs.io/en/stable/index.html

Comments

[melo_yuya]

ご回答ありがとうございます。
無事座標を抽出することが出来ました。

Answer 2

以下で取り出せると思います。

import re

def convert_coordinates(x):
    s = re.sub('(^[^(]*(*)|()*)', '', x.wkt)
    return [[float(x) for x in y.strip().split(' ')] for y in s.split(',')]

coordinates = convert_coordinates(yoshi_points[0])
print(coordinates)

Comments

[melo_yuya]

ご回答ありがとうございます。
座標を取り出すことが出来ました。
取り出した座標のうち、X座標・Y座標のみを取り出す方法があれば教えていただきたいです。

[etherbeg]

return [[float(x) for x in y.strip().split(' ')] for y in s.split(',')]
とあるところを
return [[float(x) for x in y.strip().split(' ')][:2] for y in s.split(',')]
に変えればX座標・Y座標のみを取り出すことができます。
私が回答した方法の場合は
list(p.exterior.coords)
ではなく
[[x,y] for x,y,z in p.exterior.coords]
とすればX座標・Y座標のみを取り出すことができます。

[melo_yuya]

ご回答ありがとうございます。
X座標・Y座標のみ取り出せました。