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

Question

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

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

on_road = []

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

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

Answer 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)

Answer 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