代码记录——tif网格到站点

import pandas as pd
import rasterio
import numpy as np
from scipy.ndimage import map_coordinates

# 设置Excel文件路径
stations_info_path = r'E:\2000-2018excel\黄河流域2000-2020.xlsx'

# 读取Excel文件
df = pd.read_excel(stations_info_path)

# 假设Excel文件中包含“Longitude”和“Latitude”列
lons = df['经度'].values
lats = df['纬度'].values

# 打开栅格文件
raster_path = r'E:\GSMap\daily-tif\20010720.dat.tif'
with rasterio.open(raster_path) as src:
    # 读取第一个波段的数据
    data = src.read(1)

    # 获取地理变换信息
    transform = src.transform

    # 初始化一个数组来存储插值结果
    interpolated_values = np.zeros(lons.shape)

    # 对每个站点进行插值
    for i, (lon, lat) in enumerate(zip(lons, lats)):
        # 将经纬度坐标转换为栅格坐标
        x = (lon - transform.xoff) / transform.a
        y = (transform.yoff - lat) / transform.e  # 注意：y坐标可能需要反转

        # 确保坐标在栅格范围内
        x = np.clip(x, 0, data.shape[1] - 1)
        y = np.clip(y, 0, data.shape[0] - 1)

        # 计算插值点周围的四个栅格坐标（双线性插值需要四个点）
        x0, x1 = int(np.floor(x)), int(np.ceil(x))
        y0, y1 = int(np.floor(y)), int(np.ceil(y))

        # 确保坐标不越界
        x0, x1 = np.clip([x0, x1], 0, data.shape[1] - 1)
        y0, y1 = np.clip([y0, y1], 0, data.shape[0] - 1)

        # 计算插值权重（这里简化处理，实际应考虑距离）
        fx = x - x0
        fy = y - y0

        # 双线性插值（这里使用简化的权重，实际应使用更精确的计算）
        q11 = data[y0, x0]
        q12 = data[y0, x1]
        q21 = data[y1, x0]
        q22 = data[y1, x1]

        # 注意：这里只是演示了双线性插值的基本思想，实际上你可以直接使用map_coordinates
        # 但由于map_coordinates已经内置了双线性插值，我们可以直接用它来避免手动计算

        # 使用map_coordinates进行插值（更推荐的方式）
        coords = np.array([[y], [x]])  # 注意map_coordinates需要二维数组作为输入
        interpolated_value = map_coordinates(data, coords, order=1, mode='constant', cval=np.nan).ravel()[0]

        # 存储插值结果
        interpolated_values[i] = interpolated_value

    # 现在interpolated_values包含了每个站点的插值结果
print(interpolated_values)

想进行的是tif网格数据插值到气象站点利用双线性插值，但是由于我不会写代码，很多都是搬运，双线性插值方法很复杂，老是报错，插值不成功，于是找到了这个方法。这个方法用的是线性插值，虽然不是双线性插值，但是可以作为近似。仅供参考！