パンダの緯度/経度のポイントを変換し、境界ポリゴンに入るかどうかを確認するにはどうすればいいですか？

Question

私はこのようなパンダのデータフレームdfを持っている：

id lat lon 
jhg 2.7 3.5 
ytr 3.1 3.5 
... 

私はまた、いくつかのポリゴンとGeopandasのデータフレームpolyを持っています。今、dfの点をの内側に入れて、多角形の点をプロットしたいと思います。だから、私はpoly.intersects(p)のような何かをすることができるはずです。pはかっこPointです。しかし、私は間違ったことをしています。緯度/経度のポイントでデータフレームを渡し、彼らがpolyにオーバーレイプロットするための最良の方法だろう何

from shapely.geometry import Point 
for index, row in df.iterrows(): 
    t = poly.intersects(Point(row.lon, row.lat)) 

？最小/最大緯度/経度の範囲を定義することができますが、ドットはpolyの外側にありますが、（大きな）境界ボックスの内側にも印刷されます。

Answer 1

あなたの出発点：

import pandas as pd 
from shapely.geometry import box 
import matplotlib.pyplot as plt 

from matplotlib.collections import PatchCollection 
from matplotlib.patches import Polygon 
from shapely.geometry import Point 
import seaborn as sns 
import numpy as np 

# some pretend data 
data = {'lat':[2.7,3.5,1.4,2.3,.9,1.9], 'lon':[1.2,.9,1.9,2.2,3,1.1]} 
df = pd.DataFrame(data) 

# the 'bounding' polygon 
poly = box(1,1,2,2) 
patches = PatchCollection([Polygon(poly.exterior)], facecolor='red', linewidth=.5, alpha=.5) 


# plot the bounding box 
fig, ax = sns.plt.subplots(1, figsize=(4,4)) 
ax.add_collection(patches, autolim=True) 

# plot the lat/lon points 
df.plot(x='lat',y='lon', kind='scatter',ax=ax) 
plt.show() 

の数字は次のようになります。不要なポイントを取り除くために

一つの方法は、boolean型マスクを使用することです：

#probably more efficient ways to do this, but this works 
mask = [poly.intersects(Point(lat,lon)) for lat,lon in zip(df.lat,df.lon)] 
df = df[mask] 

# make new plot (must make a new 'patch' object) 
patches1 = PatchCollection([Polygon(poly.exterior)], facecolor='red', linewidth=.5, alpha=.5) 
fig1, ax1 = sns.plt.subplots(1, figsize=(4,4)) 
ax1.add_collection(patches1, autolim=True) 

# make the axis bounds the same 
ax1.set_xlim(ax.get_xlim()) 
ax1.set_ylim(ax.get_ylim()) 

# plot the lat/lon points 
df.plot(x='lat',y='lon', kind='scatter',ax=ax1) 
plt.show() 

この画像を提供します。

あなたは緯度がポリゴンで最高点を超えているかどうかなど、他の、より高速な方法でブールマスクを作ることができることに注意してください。それらは自分自身では完璧ではないかもしれませんが、問題を減らすことができますので、intersects()を何度も呼び出す必要はありません。

[編集：ポリゴンが四角形の場合]別の方法（あなたの質問で示唆したように）は、境界ポリゴンの周りの画像を単に「トリミングする」ことです。これは、あなたがそのintersects()関数を繰り返し呼び出す必要がないので、はるかに高速なソリューションです。境界ポリゴンに基づいて画像をトリミングするには、plt.plot()前にこの権利を挿入することができます。

ax.set_xlim((np.min(poly.exterior.xy[0]),np.max(poly.exterior.xy[0]))) 
ax.set_ylim((np.min(poly.exterior.xy[1]),np.max(poly.exterior.xy[1]))) 

は、次のように与える：

Answer 2

This tutorialは、あなたがやりたいようです。また、ジオパンダの組み込みのRtree空間索引付けを利用して、高速な交差問合せを実行します。

spatial_index = gdf.sindex 
possible_matches_index = list(spatial_index.intersection(polygon.bounds)) 
possible_matches = gdf.iloc[possible_matches_index] 
precise_matches = possible_matches[possible_matches.intersects(polygon)] 

そして、あなたがやりたいように、それは、異なる色でポリゴンとその内側とその外側のポイントをプロットします。