使用 shapely 和 geopandas 检查点列表是不是接近 python 中的线串列表的快速方法

Posted 2023-03-23

【中文标题】使用 shapely 和 geopandas 检查点列表是不是接近 python 中的线串列表的快速方法

【英文标题】：Fast way to check if a list of points is close to a list of linestrings in python with shapely and geopandas使用 shapely 和 geopandas 检查点列表是否接近 python 中的线串列表的快速方法

【发布时间】：2022-01-22 22:38:44

【问题描述】：

我有大量匀称点（约 150k）和大量匀称线串（约 240k）。我想知道是否有一种快速的方法来检查这些点是否靠近任何线串。

这是我用了 300 分的代码，用时 387 秒。

on_road = []

for point in noisePoints.geometry:
    if (roads_df.geometry.distance(point) < 1e-4).any():
        on_road.append(point)

这个方法太慢了，所以我希望有一个更快的方法。

【问题讨论】：

你能分享一些数据来测试另一种方法吗？但是，我相信这实际上是一个非常简洁的解决方案。如果缓冲线串然后进行连接，则另一种选择。但我怀疑它会更慢。但一定要分享。

【参考方案1】：

您可以使用Dask-Geopandas，它可以对您的 Geopandas 数据框进行分区并并行运行任何此类匀称的操作。尝试将您的“roads_df”创建为 Dask Geopandas 数据框并尝试运行相同的代码。它应该看起来像这样

roads_ddf = dask_geopandas.from_geopandas(roads_df, npartitions=4)
on_road = []
for point in noisePoints.geometry:
    if (roads_ddf.geometry.distance(point) < 1e-4).any():
        on_road.append(point)

【参考方案2】：

一种行之有效的方法是https://geopandas.org/en/stable/docs/reference/api/geopandas.sjoin_nearest.html LineString 使用英国高速公路网络，Point 几何图形使用 NHS 医院。这提供了 227 家医院和 3669 条路段

153 ms ± 3.22 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

已使用来自条例调查的数据的 CRS。它确实适用于 epsg:4326，但是在进行空间连接时会收到警告

解决方案

# find nearest line to each point (hospital)
def nr_dist():
    gdfn = gpd.sjoin_nearest(
        gdfh.to_crs(gdf.crs),
        gdf
    ).merge(gdf["geometry"], left_on="index_right", right_index=True, suffixes=("","_line"))
    # calc the distance between the point and linestring
    gdfn["distance"] = gdfn.apply(lambda r: r["geometry_line"].distance(r["geometry"]), axis=1)
    return gdfn

# %timeit gdfn = nr_dist()
gdfn = nr_dist()

数据来源

import geopandas as gpd
import shapely.geometry
import numpy as np
import requests, io
from pathlib import Path
from zipfile import ZipFile
import pandas as pd

# get uk motorway geometry
url = "https://api.os.uk/downloads/v1/products/OpenRoads/downloads?area=GB&format=ESRI®+Shapefile&redirect"
f = Path.cwd().joinpath("road_gp.zip")

if not f.exists():
    r = requests.get(url, stream=True, headers="User-Agent": "XY")
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
    zfile = ZipFile(f)
    zfile.extractall(f.stem)

nf = Path.cwd().joinpath("a-roads")
if not nf.exists():
    nf.mkdir()
nf = nf.joinpath("a-roads.shp")
if nf.exists():
    print("exists")
    gdf = gpd.read_file(nf)
else:
    gdf = pd.concat(
        [
            gpd.read_file(sf)
            .loc[lambda d: d["class"].isin(["Motorway"])]
            .assign(source=sf.stem)
            for sf in list(Path.cwd().joinpath(f.stem).glob("**/*Link.shp"))
        ]
    )

    gdf.to_file(nf)
    
gdf = gdf.loc[~gdf["formOfWay"].isin(["Slip Road", "Roundabout"]),  ["roadNumber", "class", "source", "geometry"]]

# get hospitals and turn into geodataframe
dfhos = pd.read_csv(io.StringIO(requests.get("https://assets.nhs.uk/data/foi/Hospital.csv").text),sep="Č",engine="python",)
dfhos = dfhos.loc[lambda d: d["Sector"].eq("NHS Sector") & d["SubType"].eq("Hospital")].groupby("ParentODSCode", as_index=False).first()

# points for hospitals
gdfh = gpd.GeoDataFrame(
    geometry=dfhos.loc[:, ["Longitude", "Latitude"]].apply(
        shapely.geometry.Point, axis=1
    ), data=dfhos
).set_crs("epsg:4326").loc[:, ["Sector", "OrganisationName", "County", "geometry"]]

可视化

# visualise it...
mask = gdf.index.isin(gdfn["index_right"])
m = gdf.loc[mask].explore(color="blue", style_kwds="weight":6)
m = gdf.loc[~mask].explore(m=m, color="skyblue", style_kwds="weight":4)
m = gdfn.explore(m=m, color="red")
m