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)
  • 48
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
Python中的普通克是一种用于空间数据的方法。在这种方法中,我们使用已知的样本点数据来推断未知位置的数据。下面是一个使用Python进行普通克的示例代码: ```python import gma import pandas as pd # 读取Excel中的样本数据 Data = pd.read_excel("Interpolate.xlsx") Points = Data.loc[:, ['经度','纬度']].values Values = Data.loc[:, ['']].values # 定义参数和模型 KD = gma.smc.Interpolate.Kriging(Points, Values, Resolution=0.05, VariogramModel='Spherical', VariogramParameters=None, KMethod='Ordinary', InProjection='EPSG:4326') # 将结果写入文件 gma.rasp.WriteRaster(r'.\gma_OKriging.tif', KD.Data, Projection='WGS84', Transform=KD.Transform, DataType='Float32') ``` 以上代码中,我们首先读取了保存样本数据的Excel文件,并提取出经度和纬度的坐标以及对应的样本。然后,我们使用这些样本数据进行克,指定了的分辨率、变异函数模型等参数。最后,将结果写入一个文件中。 这样,我们就可以使用Python进行普通克了。<span class="em">1</span><span class="em">2</span><span class="em">3</span> #### 引用[.reference_title] - *1* [python:克](https://blog.csdn.net/qq_35591253/article/details/128773189)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_2"}}] [.reference_item style="max-width: 50%"] - *2* *3* [【Python进阶】克法的实现过程](https://blog.csdn.net/qq_38140292/article/details/127953822)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_2"}}] [.reference_item style="max-width: 50%"] [ .reference_list ]

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

海绵波波107

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值