Python栅格数据克里金插值

结果

输入文件

1、需要有经纬度信息以及对应的单点值

2、还要用到一个研究区的矢量文件，当然上面点的经纬度信息要在该矢量文件以内

核心代码

file_path = workspace1
    # Attempt to read the Excel file
    df = readDataFile(file_path)
    date = str("Value")
    df_cleaned = df.dropna(subset=[date])
    shp_file_path = workspace2
    region_shp = gpd.read_file(shp_file_path)
    bounds = region_shp.bounds
    data = pd.read_excel(file_path)
    coordinates = data[['POINT_X', 'POINT_Y']].values
    values = data[date].values
    longitude = data['POINT_X'].values
    latitude = data['POINT_Y'].values
    OK = OrdinaryKriging(longitude, latitude, values, variogram_model='spherical', verbose=False, enable_plotting=False)
    grid_x, grid_y = np.mgrid[bounds.minx.min():bounds.maxx.max():400j, bounds.miny.min():bounds.maxy.max():400j]
    grid_z, ss = OK.execute('grid', grid_x[:, 0], grid_y[0, :])
    transform = rasterio.transform.from_origin(grid_x.min(), grid_y.max(), (grid_x.max() - grid_x.min()) / 400,
                                               (grid_y.max() - grid_y.min()) / 400)
    out_shape = (grid_z.shape[0], grid_z.shape[1])
    mask = geometry_mask(region_shp.geometry, invert=True, transform=transform, out_shape=out_shape)
    grid_z_masked = np.where(mask, grid_z, np.nan)
    transform = from_origin(grid_x.min(), grid_y.max(), (grid_x.max() - grid_x.min()) / 400,
                            (grid_y.max() - grid_y.min()) / 400)
    with rasterio.open(workspace3, 'w', driver='GTiff',
                       height=grid_z.shape[0], width=grid_z.shape[1],
                       count=1, dtype=str(grid_z.dtype),
                       crs='+proj=latlong', transform=transform) as dst:
        dst.write(grid_z_masked, 1)
    del dst

    filename = workspace3
    with rasterio.open(filename) as src:
        # 读取第一个波段的数据
        band_data1 = src.read(1)
        # 归一化处理
        _normalized = np.floor(255 * (band_data1 - np.nanmin(band_data1)) / (np.nanmax(band_data1) - np.nanmin(band_data1)))
        # 应用色彩查找表（LUT）
        lut = np.zeros((256, 3), dtype=np.uint8)
        # [填充LUT，与原代码相同的逻辑]

        for i in range(256):
            if i < 51:
                # 红色到黄色（增加绿色）
                lut[i, 0] = i * 5
                lut[i, 1] = 150+i * 2

            elif i < 102:
                # 黄色到绿色（减少红色）
                lut[i, 0] = 255 - (i - 51) * 5
                lut[i, 1] = 255
            elif i < 153:
                # 绿色到青色（增加蓝色）
                lut[i, 1] = 255
                lut[i, 2] = (i - 102) * 5
            elif i < 204:
                # 青色到蓝色（减少绿色）
                lut[i, 1] = 255 - (i - 153) * 5
                lut[i, 2] = 255
            else:
                # 蓝色到紫色（增加红色）
                lut[i, 0] = (i - 204) * 5
                lut[i, 2] = 255
        # 应用LUT生成伪彩色图像
        colored_image = lut[np.int16(_normalized)]
        for i in range(3):
            colored_image[:, :, i] = np.where(_normalized > 0, colored_image[:, :, i], 0)

        # 创建一个新的TIFF文件来保存伪彩色图像
        with rasterio.open(
                workspace4,
                'w',
                driver='GTiff',
                height=colored_image.shape[0],
                width=colored_image.shape[1],
                count=3,
                dtype=colored_image.dtype,
                crs=src.crs,
                transform=src.transform,
        ) as dst:
            for i in range(3):
                dst.write(colored_image[:, :, i], i + 1)