克里金krige

什么是克里金（krige）

插值模式

• 普通克里金插值
• 泛克里金插值

半变异模型

• 球面函数
• 三角函数
• 指数函数
• 高斯函数
• 线性函数

输出图像元大小

当前进度

由于插值数量太大，经常导致内存不足，无法继续

import numpy as np
from pykrige.ok import OrdinaryKriging
from shapely.geometry import Point
from shapely.ops import unary_union
import geopandas as gpd


import geopandas as gpd
import numpy as np
import pandas as pd

from pykrige.ok import OrdinaryKriging
from scipy.spatial import KDTree


def kriging(
    input_shp,
    data_field,
    resolution: int = 100,  # 采用百分比，100则为100分之一，1000则为1000份之一
    semivariogram_model="spherical",
    n_nearest_points=12,
    range_max=None,
):
    # 读取输入数据
    gdf = gpd.read_file(input_shp)
    values = gdf[data_field].values
    # 使用np.meshgrid生成网格的X和Y坐标
    X, Y = np.meshgrid(grid_x, grid_y)

    # 将X和Y坐标展平为一维数组，因为PyKrige需要这样的格式
    X_flat = X.flatten()
    Y_flat = Y.flatten()

    # 初始化克里金插值器
    OK = OrdinaryKriging(
        coordinates[:, 0],
        coordinates[:, 1],
        values,
        variogram_model=semivariogram_model,
        verbose=False,
        enable_plotting=False,
        nlags=n_nearest_points,
        range_max=range_max,
    )

    # 执行克里金插值
    z, ss = OK.execute("grid", X_flat, Y_flat)

    # 将结果重新整形为二维数组，以匹配网格的形状
    z_reshaped = z.reshape(X.shape)

    # 返回插值结果
    return z_reshaped


def extract_kriging_results(kriging_result, points_shp, output_field="kriging_value"):
    # 读取点集shp文件
    points_gdf = gpd.read_file(points_shp)

    # 创建一个空列表来存储每个点的插值结果
    kriging_values = []

    # 遍历点集中的每个点
    for index, row in points_gdf.iterrows():
        x, y = row.geometry.x, row.geometry.y

        # 计算该点在kriging结果中的索引（这里假设kriging_result是一个二维numpy数组）
        # 注意：这里使用简单的线性插值索引计算，可能需要根据实际情况调整
        # （例如，如果kriging_result的分辨率不是严格的output_resolution，则需要进行适当的缩放）
        ix = int((x - kriging_result.shape[1] * output_resolution / 2) / output_resolution)
        iy = int((y - kriging_result.shape[0] * output_resolution / 2) / output_resolution)

        # 检查索引是否在范围内
        if 0 <= ix < kriging_result.shape[1] and 0 <= iy < kriging_result.shape[0]:
            # 提取插值结果
            kriging_values.append(kriging_result[iy, ix])
        else:
            # 如果点不在kriging网格内，可以选择一个默认值（例如NaN）
            kriging_values.append(np.nan)

    # 将插值结果作为新列添加到GeoDataFrame中
    points_gdf[output_field] = kriging_values

    # 输出新的shp点文件
    output_shp = "output_points_with_kriging.shp"
    points_gdf.to_file(output_shp, driver="ESRI Shapefile")

    # 返回更新后的GeoDataFrame（可选）
    return points_gdf


if __name__ == "__main__":
    input_shp = r"./clip_range.shp"
    point_shp = r"./fishnet_clip.shp"
    af_shp = r"./af_raw.shp"
    be_shp = r"./be_raw.shp"

    gdf = gpd.read_file(be_shp)
    Z = gdf["speed"].values
    xmin, ymin, xmax, ymax = gdf.total_bounds

    x_l = np.linspace(xmin, xmax, 300)
    y_l = np.linspace(ymin, ymax, 100)

    # 将X和Y坐标展平为一维数组，因为PyKrige需要这样的格式
    xgrid, ygrid = np.meshgrid(x_l, y_l)

    coordinates = np.array([(row.geometry.x, row.geometry.y) for _, row in gdf.iterrows()])
    X = gdf["geometry"].apply(lambda geom: geom.x).values
    Y = gdf["geometry"].apply(lambda geom: geom.y).values
    OK = OrdinaryKriging(X, Y, Z, variogram_model="gaussian", nlags=12)
    zgrid, ss = OK.execute("grid", x_l, y_l)

    print(len(zgrid))
    print(len(ygrid))