确定一个点在哪个多边形中，然后将该多边形的名称作为一个新列应用到一个大熊猫数据框

Posted 2023-04-18

【中文标题】确定一个点在哪个多边形中，然后将该多边形的名称作为一个新列应用到一个大熊猫数据框

【英文标题】：Determine in which polygon a point is, and then apply the name of that polygon as a new column to a large pandas dataframe

【发布时间】：2020-03-26 17:04:25

【问题描述】：

我有一个大型数据框，其中包含来自世界各地各种船只的位置数据。 imoNo是船舶标识符。以下是数据框的示例： 这是重现它的代码：

# intialise data of lists. 
ships_data = 'imoNo':[9321483, 9321483, 9321483, 9321483, 9321483], 
                    'Timestamp':['2020-02-22 00:00:00', '2020-02-22 00:10:00', '2020-02-22 00:20:00', '2020-02-22 00:30:00', '2020-02-22 00:40:00'],
                    'Position Longitude':[127.814598, 127.805634, 127.805519, 127.808548, 127.812370],
                    'Position Latitude':[33.800232, 33.801899, 33.798885, 33.795799, 33.792931] 

# Create DataFrame 
ships_df = pd.DataFrame(ships_data) 

我需要做的是在数据框的末尾添加一列，该列将标识船舶航行的sea_name。类似于以下内容：

为了到达那里，我在this 链接（IHO Sea Areas v3）找到了一个 .shp 格式的数据集，如下所示：

所以，我这样做的方式是遍历船舶数据集的每个 (long,lat)，检查该对在哪个多边形中，最后返回匹配多边形的海洋名称。这是我的代码：

### Load libraries
import numpy as np
import pandas as pd
import geopandas as gp
import shapely.speedups
from shapely.geometry import Point, Polygon
shapely.speedups.enable()

### Check and map lon lat pair with sea name
def get_seaname(long,lat):
    pnt = Point(long,lat)
    for i,j in enumerate(iho_df.geometry):
        if pnt.within(j):
            return iho_df.NAME.iloc[i]


### Apply the above function to the dataframe
ships_df['sea_name'] = ships_df.apply(lambda x: get_seaname(x['Position Longitude'], x['Position Latitude']), axis=1) 

但是，这是一个非常耗时的过程。我在 Mac 上的前 1000 行 ship_df 上进行了本地测试，运行大约需要 1 分钟。如果它呈线性增长，那么整个数据集大约需要 14 天：-D。

任何优化上述功能的想法将不胜感激。

谢谢！

【问题讨论】：

有趣的问题，但是您能否将您的示例数据发布为文本而不是图像？使它更加更容易工作。

有趣.. 我当然可以应用矢量化来加快计算速度，但我想知道如何加快将 lat 和 long 转换为 namedtuple（我的矢量化想法的先决条件）而不让它变得愚蠢慢... :(

【参考方案1】：

最后，我比最初的问题更快。

首先，我使用来自 IHO 海域数据集的信息创建了一个描述边界框的多边形

# Create a bbox polygon
iho_df['bbox'] = iho_df.apply(lambda x: Polygon([(x['min_X'], x['min_Y']), (x['min_X'], x['max_Y']), (x['max_X'], x['max_Y']), (x['max_X'], x['min_Y'])]), axis=1)

然后，我更改了函数，以便首先查看 bbox（这比几何图形快得多，因为它只是一个矩形）。当一个点落在多个框内（用于接壤的海洋）时，它才会查看初始多边形以在匹配框（而不是所有多边形内）中找到正确的海洋名称。

# Function that checks and maps lon lat pair with sea name
def get_seaname(long,lat):
    pnt = Point(long,lat)
    names = []
    # Check within each bbox first to note the polygons to look at
    for i,j in enumerate(iho_df.bbox):
        if pnt.within(j):
            names.append(iho_df.NAME.iloc[i])
    # Return nan for no return
    if len(names)==0:
        return np.nan
    # Return the single name of the sea 
    elif len(names)==1:
        return names[0]
    # Run the pnt.within() only for the polygons within the collected sea names
    else:
        limitizez_df = iho_df[iho_df['NAME'].isin(names)].reset_index(drop=True)
        for k,n in enumerate(limitizez_df.geometry):
            if pnt.within(n):
                return limitizez_df.NAME.iloc[k]

这个大大减少了时间。为了进一步提升它，我使用多处理在 CPU 内核之间进行并行化。这个想法来自另一个 *** 帖子，我现在不记得了，但这里是代码。

import multiprocessing as mp

# Function that parallelizes the apply function among the cores of the CPU
def parallelize_dataframe(df, func, n_cores):
    df_split = np.array_split(df, n_cores)
    pool = Pool(n_cores)
    df = pd.concat(pool.map(func, df_split))
    pool.close()
    pool.join()
    return df

# Function that adds a sea_name column in the main dataframe
def add_features(df):
    # Apply the function
    df['sea_name'] = df.apply(lambda x: get_seaname(x['Position Longitude'], x['Position Latitude']), axis=1)
    return df

最后，我没有对 get_seaname() 使用 apply 函数，而是将它用于 parallelize_dataframe() 函数，以在所有可用的 CPU 内核上运行：

### Apply the above function to the dataframe
ships_df = parallelize_dataframe(ships_df, add_features, n_cores=mp.cpu_count())

我希望我的解决方案也能帮助其他人！

【参考方案2】：

试试这个，

用apply 找出每个lat long （无法实现更快的方法，欢迎帮助）

import numpy as np
import pandas as pd
import geopandas as gp
import shapely.speedups
from shapely.geometry import Point, Polygon
shapely.speedups.enable()

# I am still uncomfortable with this. More ideas on speeding up this part are welcome
ships_df['point'] = ships_df.apply(lambda x: Point(x['Position Longitude'], x['Position Latitude']), axis=1)

现在矢量化您的函数以在 Point 上工作

def get_seaname(pnt:Point):
    for i,j in enumerate(iho_df.geometry):
        if pnt.within(j):
            return iho_df.NAME.iloc[i]

现在，由于您的方法适用于单个点，请将点列转换为 Point 对象的向量并向量化您的方法

get_seaname = np.vectorize(get_seaname)

ships_df['sea_name'] = pd.Series(get_seaname(ships_df['point'].values))

【讨论】：

感谢您的快速回复！我在 1000 行样本上尝试了你的建议，但仍然需要很长时间（大约 1 分钟），在完成之前我收到一个错误：/Users/oikonang/miniconda3/envs/geo/lib/python3.7/site-packages/numpy/lib/function_base.py:2167: RuntimeWarning: invalid value encountered in ? (vectorized) outputs = ufunc(*inputs)