Geopandas 检查点是否在多边形内

Question

我有海洋 geopandas，其中包含 1 个多边形（来源：naturalearthdata.com）

我还有另一个包含大量经度和纬度信息的数据框

我想添加一个新列，如果该点位于海洋中（多边形内部），则该列将为 True

zipfile = "ne_10m_ocean/ne_10m_ocean.shp"
ocean_gpd = geopandas.read_file(zipfile)

df = pd.DataFrame({
    'lon': [120.0,120.1,120.2,120.3,120.4],
    'lat': [10.0,10.1,10.2,10.3,10.4]
})

for index, row in df.iterrows():
    df.loc[index,'is_ocean'] = ocean_gpd.contains(Point(x['lon'],x['lat'])

但是太慢了，我尝试使用这样的 lambda 函数

df = df.assign(is_ocean = lambda x: ocean_gpd.contains(Point(x['lon'],x['lat']))

但是失败了，错误是

cannot convert the series to <class 'float'>

有人知道如何在 geopandas 中进行更好的个人点检查吗？

注意：我刚刚意识到对于多边形数据我使用了10m一个（更详细的多边形），如果我使用110m会好很多，但将来也许我需要使用10m

Answer 1

您可以像这样使用

apply

：

import geopandas
import pandas as pd
from shapely.geometry import Point

ocean_gpd = geopandas.read_file('ne_10m_ocean.shp')

df = pd.DataFrame({
    'lon': [120.0, 120.1, 120.2, 120.3, 120.4],
    'lat': [10.0, 10.1, 10.2, 10.3, 10.4]
})

def in_ocean(row):
    point = Point(row['lon'], row['lat'])
    return ocean_gpd.contains(point).any()

df['is_ocean'] = df.apply(in_ocean, axis=1)

返回：

     lon   lat  is_ocean
0  120.0  10.0     False
1  120.1  10.1     False
2  120.2  10.2     False
3  120.3  10.3     False
4  120.4  10.4     False

Answer 2

接受的答案都是上帝，尽管有很多点，它变得非常慢。

两种加速可能性是

• 并行化
• 将大多边形切碎（即一个海洋多边形切成较小的多边形）

我只用 100 分进行测试，标准方法的速度明显减慢。对整个地图上的随机点进行采样并运行代码。

import geopandas
import pandas as pd
from shapely.geometry import Point
import random
ocean_gpd = geopandas.read_file('~/ne/ne_10m_ocean.shp')
npoints = 100
random.seed(42)
df = pd.DataFrame({
    'lon': [(random.random()*360)-180 for i in range(npoints)],
    'lat': [(random.random()*180)-90 for i in range(npoints)],
})

现在

df['is_ocean'] = df.apply(in_ocean, axis=1)

在我的旧但仍然强大的 32 核服务器上大约需要 1 分钟。

如果我们更多地利用

geopandas

并创建我们自己的并行化方法，我们至少可以获得 2 倍的加速。

首先我们可以用点创建一个合适的 GeoDataFrame

gdf = geopandas.GeoDataFrame(df[['lat','lon']],
          geometry=geopandas.points_from_xy(df['lon'],df['lat']), 
          crs=ocean_gpd.crs
      )

接下来我们也可以使用

geopandas

工具来绘制数据。

a = ocean_gpd.plot()
gdf.plot(ax=a, color='red')

接下来我们需要一个函数来并行化事物。这是我并行化尚未优化的事物的首选解决方案。

from joblib import parallel_backend, Parallel, delayed, effective_n_jobs
from sklearn.utils import gen_even_slices
from sklearn.utils.validation import _num_samples

def parallel_apply(df, func, n_jobs= -1, **kwargs):
    """ Pandas apply in parallel using joblib. 
    Uses sklearn.utils to partition input evenly.
    
    Args:
        df: Pandas DataFrame, Series, or any other object that supports slicing and apply.
        func: Callable to apply
        n_jobs: Desired number of workers. Default value -1 means use all available cores.
        **kwargs: Any additional parameters will be supplied to the apply function
        
    Returns:
        Same as for normal Pandas DataFrame.apply()
        
    """
    
    if effective_n_jobs(n_jobs) == 1:
        return df.apply(func, **kwargs)
    else:
        ret = Parallel(n_jobs=n_jobs)(
            delayed(type(df).apply)(df[s], func, **kwargs)
            for s in gen_even_slices(_num_samples(df), effective_n_jobs(n_jobs)))
        return pd.concat(ret)

现在我们可以并行检查：

gdf['is_ocean'] = parallel_apply(gdf, lambda x: ocean_gpd.contains(x.geometry).any(), axis=1)

对我来说现在大约需要 30 秒。切碎大多边形将使速度更快，因为它可以使用空间树。以前我使用过此页面中的“武士刀”https://snorfalorpagus.net/blog/2016/03/13/splitting-large-polygons-for-faster-intersections/

请注意，上面的并行化函数有时似乎并没有并行化事物。我不知道这是否与 GeoPandas 或其他有关。