150 lines
4.1 KiB
Python
150 lines
4.1 KiB
Python
"""
|
||
geo_tools.analysis.spatial_ops
|
||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
空间叠加与邻域分析操作。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
from typing import Any
|
||
|
||
import geopandas as gpd
|
||
import pandas as pd
|
||
|
||
from app.utils.logger import get_logger
|
||
|
||
logger = get_logger(__name__)
|
||
|
||
|
||
def buffer_and_overlay(
|
||
source: gpd.GeoDataFrame,
|
||
distance: float,
|
||
target: gpd.GeoDataFrame,
|
||
how: str = "intersection",
|
||
projected_crs: str | None = None,
|
||
) -> gpd.GeoDataFrame:
|
||
"""对 source 执行缓冲区后与 target 执行叠置分析。
|
||
|
||
Parameters
|
||
----------
|
||
source:
|
||
源图层（生成缓冲区）。
|
||
distance:
|
||
缓冲距离（与 ``projected_crs`` 单位一致）。
|
||
target:
|
||
叠置目标图层。
|
||
how:
|
||
叠置类型：``"intersection"``、``"union"``、``"difference"``、``"symmetric_difference"``、``"identity"``。
|
||
projected_crs:
|
||
执行缓冲区前先投影到此 CRS（建议使用平面坐标系以保证距离精度）；
|
||
``None`` 则使用 source 的当前 CRS（地理 CRS 下 distance 单位为度）。
|
||
|
||
Returns
|
||
-------
|
||
gpd.GeoDataFrame
|
||
"""
|
||
original_crs = source.crs
|
||
|
||
if projected_crs:
|
||
source = source.to_crs(projected_crs)
|
||
target = target.to_crs(projected_crs)
|
||
|
||
buffered = source.copy()
|
||
buffered["geometry"] = buffered.geometry.buffer(distance)
|
||
logger.debug("缓冲区完成（distance=%.2f），执行叠置分析（how=%s）", distance, how)
|
||
|
||
result = gpd.overlay(buffered, target, how=how, keep_geom_type=False)
|
||
|
||
if projected_crs:
|
||
result = result.to_crs(original_crs) # type: ignore
|
||
|
||
logger.info("叠置分析完成：%d 条结果", len(result))
|
||
return result
|
||
|
||
|
||
def overlay(
|
||
df1: gpd.GeoDataFrame,
|
||
df2: gpd.GeoDataFrame,
|
||
how: str = "intersection",
|
||
keep_geom_type: bool = True,
|
||
) -> gpd.GeoDataFrame:
|
||
"""封装 geopandas overlay，自动对齐 CRS。
|
||
|
||
Parameters
|
||
----------
|
||
how:
|
||
叠置类型：``"intersection"``、``"union"``、``"difference"``、
|
||
``"symmetric_difference"``、``"identity"``。
|
||
"""
|
||
if df1.crs != df2.crs:
|
||
df2 = df2.to_crs(df1.crs) # type: ignore
|
||
result = gpd.overlay(df1, df2, how=how, keep_geom_type=keep_geom_type)
|
||
logger.debug("overlay(%s)：%d 条结果", how, len(result))
|
||
return result
|
||
|
||
|
||
def nearest_features(
|
||
source: gpd.GeoDataFrame,
|
||
target: gpd.GeoDataFrame,
|
||
k: int = 1,
|
||
max_distance: float | None = None,
|
||
) -> gpd.GeoDataFrame:
|
||
"""为 source 中每条要素找到 target 中最近的 k 个要素。
|
||
|
||
Parameters
|
||
----------
|
||
source:
|
||
查询图层。
|
||
target:
|
||
被查询图层。
|
||
k:
|
||
最近邻数量。
|
||
max_distance:
|
||
最大搜索距离（与 CRS 单位一致），``None`` 表示无限制。
|
||
|
||
Returns
|
||
-------
|
||
gpd.GeoDataFrame
|
||
连接了最近 target 属性的 source GDF（可能包含重复行，每行对应一个近邻）。
|
||
"""
|
||
if source.crs != target.crs:
|
||
target = target.to_crs(source.crs) # type: ignore
|
||
|
||
result = gpd.sjoin_nearest(
|
||
source,
|
||
target,
|
||
how="left",
|
||
max_distance=max_distance,
|
||
distance_col="nearest_distance",
|
||
lsuffix="left",
|
||
rsuffix="right",
|
||
)
|
||
logger.debug("最近邻分析完成（k=%d）：%d 条结果", k, len(result))
|
||
return result
|
||
|
||
|
||
def select_by_location(
|
||
source: gpd.GeoDataFrame,
|
||
selector: gpd.GeoDataFrame,
|
||
predicate: str = "intersects",
|
||
) -> gpd.GeoDataFrame:
|
||
"""按位置关系从 source 中选取要素（等同于 ArcGIS「按位置选择」）。
|
||
|
||
Parameters
|
||
----------
|
||
predicate:
|
||
空间谓词：``"intersects"``、``"within"``、``"contains"``、``"touches"``。
|
||
|
||
Returns
|
||
-------
|
||
gpd.GeoDataFrame
|
||
满足条件的 source 子集。
|
||
"""
|
||
if source.crs != selector.crs:
|
||
selector = selector.to_crs(source.crs) # type: ignore
|
||
|
||
joined = gpd.sjoin(source, selector, how="inner", predicate=predicate)
|
||
result = source.loc[source.index.isin(joined.index)].copy()
|
||
logger.debug("按位置选择（%s）：%d / %d 条", predicate, len(result), len(source))
|
||
return result
|