@step(id="analysis.voronoi", ...)
def voronoi(ctx):
    # 延迟导入——仅在步骤执行时才检查依赖
    from scipy.spatial import Voronoi as ScipyVoronoi
    from shapely.geometry import Polygon
    ...

from scipy.spatial import Voronoi

# 输入点坐标
points = np.array([
    [0, 0], [1, 0], [0, 1], [1, 1], [0.5, 0.5]
])

# 计算 Voronoi 图
vor = Voronoi(points)

# vor.vertices: Voronoi 顶点坐标
# vor.regions: 每个区域的顶点索引列表
# vor.ridge_vertices: 每条边的顶点索引对
# vor.point_region: 每个输入点对应的区域索引

# buffer=0: 裁剪到点集的精确边界框
# buffer=100: 边界框向外扩展 100 个单位

@step(id="analysis.voronoi", ...)
def voronoi(ctx):
    from scipy.spatial import Voronoi as ScipyVoronoi
    from shapely.geometry import Polygon, box
    from shapely.ops import clip_by_rect

    gdf = ctx.input("input")
    clip_to_bounds = ctx.param("clip_to_bounds")
    buffer_dist = ctx.param("buffer") or 0.0

    if clip_to_bounds is None:
        clip_to_bounds = True

    # 提取点坐标
    coords = np.array([(geom.x, geom.y) for geom in gdf.geometry])

    # 计算 Voronoi 图
    vor = ScipyVoronoi(coords)

    # 计算边界框
    minx, miny = coords.min(axis=0) - buffer_dist
    maxx, maxy = coords.max(axis=0) + buffer_dist
    bounding_box = box(minx, miny, maxx, maxy)

    # 构建多边形
    polygons = []
    for point_idx, region_idx in enumerate(vor.point_region):
        region = vor.regions[region_idx]
        if -1 in region or len(region) == 0:
            # 无界区域：使用大多边形近似后裁剪
            polygon = _construct_unbounded_region(vor, point_idx, bounding_box)
        else:
            vertices = [vor.vertices[i] for i in region]
            polygon = Polygon(vertices)

        if clip_to_bounds and polygon is not None:
            polygon = polygon.intersection(bounding_box)

        polygons.append(polygon)

    # 构建输出 GeoDataFrame
    result = gpd.GeoDataFrame(
        gdf.drop(columns="geometry").copy(),
        geometry=polygons,
        crs=gdf.crs,
    )

    stats = {
        "feature_count": len(result),
        "input_point_count": len(gdf),
    }

    return StepResult(output=result, stats=stats)

@step(id="analysis.heatmap", ...)
def heatmap(ctx):
    from scipy.stats import gaussian_kde

    gdf = ctx.input("input")
    resolution = ctx.param("resolution") or 100
    radius = ctx.param("radius")
    weight_field = ctx.param("weight_field")

    # 提取坐标
    coords = np.array([(geom.x, geom.y) for geom in gdf.geometry])
    x, y = coords[:, 0], coords[:, 1]

    # 提取权重
    weights = None
    if weight_field and weight_field in gdf.columns:
        weights = gdf[weight_field].values

    # 构建网格
    xmin, ymin, xmax, ymax = gdf.total_bounds
    xi = np.linspace(xmin, xmax, resolution)
    yi = np.linspace(ymin, ymax, resolution)
    xx, yy = np.meshgrid(xi, yi)
    grid_coords = np.vstack([xx.ravel(), yy.ravel()])

    # 核密度估计
    if radius is not None:
        bw_method = radius / np.std(coords, axis=0).mean()
    else:
        bw_method = "scott"

    kde = gaussian_kde(
        np.vstack([x, y]),
        bw_method=bw_method,
        weights=weights,
    )
    density = kde(grid_coords).reshape(resolution, resolution)

    # 翻转 y 轴（栅格的 y 轴从上到下）
    density = np.flipud(density)

    # 构建 raster_info
    pixel_width = (xmax - xmin) / resolution
    pixel_height = (ymax - ymin) / resolution
    transform = rasterio.transform.from_bounds(
        xmin, ymin, xmax, ymax, resolution, resolution
    )

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "width": resolution,
        "height": resolution,
        "count": 1,
        "crs": gdf.crs,
        "transform": transform,
        "nodata": 0.0,
    }

    result = {
        "data": density[np.newaxis, ...].astype("float32"),
        "transform": transform,
        "crs": gdf.crs,
        "profile": profile,
        "path": None,
    }

    stats = {
        "resolution": resolution,
        "density_min": float(density.min()),
        "density_max": float(density.max()),
        "point_count": len(gdf),
    }

    return StepResult(output=result, stats=stats)

from scipy.interpolate import griddata

# points: 采样点坐标 (N, 2)
# values: 采样点值 (N,)
# grid_points: 网格坐标 (M, 2)

result = griddata(
    points,
    values,
    grid_points,
    method="linear",  # 或 "nearest", "cubic"
)

# IDW 自定义实现
def idw_interpolate(points, values, grid_points, power=2.0):
    result = np.zeros(len(grid_points))

    for i, (gx, gy) in enumerate(grid_points):
        distances = np.sqrt(
            (points[:, 0] - gx)**2 + (points[:, 1] - gy)**2
        )

        # 避免除零
        zero_mask = distances == 0
        if zero_mask.any():
            result[i] = values[zero_mask][0]
        else:
            weights = 1.0 / distances ** power
            result[i] = np.sum(weights * values) / np.sum(weights)

    return result

@step(id="analysis.interpolate", ...)
def interpolate(ctx):
    from scipy.interpolate import griddata

    gdf = ctx.input("input")
    value_field = ctx.param("value_field")
    method = ctx.param("method") or "linear"
    resolution = ctx.param("resolution") or 100
    power = ctx.param("power") or 2.0

    # 提取坐标和值
    coords = np.array([(geom.x, geom.y) for geom in gdf.geometry])
    values = gdf[value_field].values.astype(float)

    # 构建网格
    xmin, ymin, xmax, ymax = gdf.total_bounds
    xi = np.linspace(xmin, xmax, resolution)
    yi = np.linspace(ymin, ymax, resolution)
    xx, yy = np.meshgrid(xi, yi)

    if method == "idw":
        # 自定义 IDW 实现
        grid_points = np.column_stack([xx.ravel(), yy.ravel()])
        result = _idw(coords, values, grid_points, power)
        grid_z = result.reshape(resolution, resolution)
    else:
        # 使用 scipy griddata
        grid_z = griddata(
            coords, values,
            (xx, yy),
            method=method,
        )

    # 翻转 y 轴
    grid_z = np.flipud(grid_z)

    # 处理 NaN（外推区域）
    grid_z = np.nan_to_num(grid_z, nan=0.0)

    # 构建 raster_info
    transform = rasterio.transform.from_bounds(
        xmin, ymin, xmax, ymax, resolution, resolution
    )

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "width": resolution,
        "height": resolution,
        "count": 1,
        "crs": gdf.crs,
        "transform": transform,
        "nodata": 0.0,
    }

    raster_info = {
        "data": grid_z[np.newaxis, ...].astype("float32"),
        "transform": transform,
        "crs": gdf.crs,
        "profile": profile,
        "path": None,
    }

    stats = {
        "method": method,
        "resolution": resolution,
        "value_min": float(np.nanmin(grid_z)),
        "value_max": float(np.nanmax(grid_z)),
        "sample_count": len(gdf),
    }

    return StepResult(output=raster_info, stats=stats)

# 聚类前
│ id │ name   │ geometry      │
├────┼────────┼───────────────┤
│ 1  │ 点A    │ POINT(1, 2)   │
│ 2  │ 点B    │ POINT(1.1, 2) │
│ 3  │ 点C    │ POINT(5, 5)   │

# DBSCAN 聚类后
│ id │ name   │ geometry      │ cluster_id │
├────┼────────┼───────────────┼────────────┤
│ 1  │ 点A    │ POINT(1, 2)   │ 0          │
│ 2  │ 点B    │ POINT(1.1, 2) │ 0          │
│ 3  │ 点C    │ POINT(5, 5)   │ -1         │  ← 噪声点

@step(id="analysis.cluster", ...)
def cluster(ctx):
    from sklearn.cluster import DBSCAN, KMeans

    gdf = ctx.input("input")
    method = ctx.param("method") or "dbscan"
    n_clusters = ctx.param("n_clusters") or 5
    eps = ctx.param("eps") or 0.5
    min_samples = ctx.param("min_samples") or 5

    # 提取坐标
    coords = np.array([(geom.x, geom.y) for geom in gdf.geometry])

    # 执行聚类
    if method == "dbscan":
        model = DBSCAN(eps=eps, min_samples=min_samples)
        labels = model.fit_predict(coords)
    elif method == "kmeans":
        model = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)
        labels = model.fit_predict(coords)
    else:
        raise ValueError(f"不支持的聚类方法: {method}")

    # 添加 cluster_id 列
    result = gdf.copy()
    result["cluster_id"] = labels

    # 统计信息
    unique_labels = set(labels)
    n_clusters_found = len(unique_labels - {-1})
    n_noise = int((labels == -1).sum())

    stats = {
        "feature_count": len(result),
        "method": method,
        "n_clusters_found": n_clusters_found,
        "n_noise_points": n_noise,
    }

    return StepResult(output=result, stats=stats)