将点捕捉到形状中最近点处的多边形/线

Question

我有一个点层（点）和一个线层（子网格），它是将多边形边界转换为线。我希望这些点捕捉到线上最近的点（如果在 100m 之内），这样我就可以使用空间连接从中获取属性。

我尝试过使用 shapely.ogr.snap 但它非常关闭，就像这样

绿点应该捕捉到右侧而不是左侧的多边形 我尝试遵循涉及插值/项目的其他答案，但我得到“AttributeError：'GeoDataFrame'对象没有属性'_geom'”。我需要使用 wkt.loads 或 unary_union 吗？到目前为止，这些对我来说都失败了。

point4 = point3.copy()point4['geometry'] = point4['geometry'].astype(str).apply(wkt.loads)
point4 = point4.set_geometry(col='geometry')

subgrid2 = geopandas.read_file(f"id_{545}.gpkg")
subgrid2['gridcell'] = subgrid2.id.astype(int)
subgrid2 = subgrid2[['gridcell', 'geometry']]
subgrid2.set_geometry('geometry', inplace=True)
subgrid2 = subgrid2.to_crs(epsg=2278)
subgrid2['geometry'] = subgrid2.geometry.astype(str).apply(wkt.loads)
#subgrid2 = subgrid2.unary_union

wkt.loads(str(subgrid2.iloc[0, 1])))[2].wkt
#geopandas.GeoDataFrame(subgrid2.interpolate(subgrid2.project(point4.geometry))).to_csv("fff.csv")

Answer 1

使用 geopandas.sjoin_nearest 查找 100 以内最近的线路米。

shapely.ops.nearest_points 查找上最近的点最近的线路。

然后 shapely.ops.snap 到它：

import geopandas as gpd
from shapely.ops import nearest_points, snap

line = gpd.read_file(r"C:\GIS\data\testdata\6k_lines.geojson")
point = gpd.read_file(r"C:\GIS\data\testdata\1000points.geojson")

#Create unique id
point["pointid"] = range(point.shape[0])

line["linegeom"] = line.geometry #Save the line geometry or it is lost in the spatial join

sj = gpd.sjoin_nearest(left_df=point, right_df=line, how="left", max_distance=100)
#print(sj.pointid.duplicated().any()) #If a point is exactly at the same distance to two lines, it can be duplicated
#False #OK

# sj.dtypes
# geometry       geometry
# pointid           int64
# index_right     float64
# lineid          float64
# linegeom       geometry

#For rows with geometries in linegeom columns, the point should be snapped to the closest point on it:
#sj.head()
#                          geometry  ...                                           linegeom
# 0  POINT (558376.685 7048618.235)  ...                                               None
# 1  POINT (534507.029 6997811.101)  ...                                               None
# 2  POINT (546192.206 7004679.200)  ...                                               None
# 3  POINT (551066.576 7007449.539)  ...  LINESTRING (550000.000 7007500.000, 552500.000...

#Find closest point
sj["closest_point"] = sj.apply(lambda x: nearest_points(x.geometry, x.linegeom)[1] if x.linegeom is not None else None, axis=1)

#If there is a point to snap to, snap. Else leave the original geometry intact.
sj["geometry"] = sj.apply(lambda x: snap(x.geometry, x.closest_point, 110) if x.closest_point is not None else x.geometry, axis=1)

sj = sj[[col for col in sj.columns if col in point.columns]] #Drop all columns created in the processing

sj.to_file(r"C:\GIS\data\testdata\1000points_SNAPPED.geojson")

最左边和最右边的点被捕捉到线条。其余的都不是因为距离> 100米。