学习记录3——利用python使用双线性插值方法将nc格点数据插值到气象站点

import os
import pandas as pd
import numpy as np
import netCDF4 as nc

# 站点信息，包含区站号、经度、纬度
stations_info = pd.read_excel(r'E:\2000-2018excel\黄河流域2000-2020.xlsx', 'Sheet1',
                              index_col=None, na_values=['NA'])

lonSta = stations_info['经度'].to_numpy()
latSta = stations_info['纬度'].to_numpy()

# 定义双线性插值函数
def Bilinear_interp(lonSta, latSta, longitude, latitude, var):
    var_sta = np.zeros((len(lonSta), 1))
    for i in range(len(lonSta)):
        # 找到四个最近邻点的索引
        iSta = np.searchsorted(longitude, lonSta[i]) - 1
        jSta = np.searchsorted(latitude, latSta[i]) - 1

        # 边界条件处理
        iSta = np.clip(iSta, 0, len(longitude) - 1)
        jSta = np.clip(jSta, 0, len(latitude) - 1)

       # print("iSta:", iSta, "jSta:", jSta)
        #print("var shape:", var.shape)
        # 获取四个插值点的值
        var11 = var[jSta, iSta]
        var21 = var[jSta, min(iSta + 1, len(longitude) - 1)]
        var12 = var[min(jSta + 1, len(latitude) - 1), iSta]
        var22 = var[min(jSta + 1, len(latitude) - 1), min(iSta + 1, len(longitude) - 1)]

        valid_values = [v for v in [var11, var21, var12, var22] if not np.isnan(v)]
        # 检查四个值是否全为NaN或零
        if all(np.isnan(v) or v == 0 for v in [var11, var21, var12, var22]):
            var_interp = 0  # 如果所有值都为NaN或零，则插值为0
        elif len(valid_values) < 4:
            # 如果少于4个有效值，使用最邻近的有效值
            nearest_value = valid_values[0]  # 这里简化为只取第一个有效值
            var_interp = nearest_value
        else:
            # 双线性插值计算
            x = lonSta[i]
            y = latSta[i]
            x1 = longitude[iSta]
            x2 = longitude[min(iSta + 1, len(longitude) - 1)]
            y1 = latitude[jSta]
            y2 = latitude[min(jSta + 1, len(latitude) - 1)]
            if (x2 - x1) == 0 or (y2 - y1) == 0:
                nearest_value = valid_values[0] if valid_values else np.nan
                var_interp = nearest_value
            else:
                arg = 1.0 / ((x2 - x1) * (y2 - y1))
                arg1 = (x2 - x) * (y2 - y) * arg
                arg2 = (x - x1) * (y2 - y) * arg
                arg3 = (x2 - x) * (y - y1) * arg
                arg4 = (x - x1) * (y - y1) * arg
                var_interp = var11 * arg1 + var21 * arg2 + var12 * arg3 + var22 * arg4

        var_sta[i] = var_interp

    return var_sta.ravel()

# 处理输入目录中的所有NetCDF文件
input_dir = r'E:\satlite\MSWEP'
for filename in os.listdir(input_dir):
    if filename.endswith('.nc'):

        input_path = os.path.join(input_dir, filename)

        # 打开NetCDF文件
        dataset = nc.Dataset(input_path)
        #读取变量的维度
        for varname in dataset.variables:
            print(varname, dataset.variables[varname].dimensions)

        # 读取经纬度和降水数据
        longitude = dataset.variables['lon'][:].data
        latitude = dataset.variables['lat'][:].data
        #print(latitude)
        #如果纬度是从大到小，直接插值会出现错误，需要将纬度翻转
        latitude = latitude[::-1]  # 翻转纬度数组
        #print(latitude)
        #二维
        #precipitation = dataset.variables['Band1'].data
        #三维
        precipitation = dataset.variables['precipitation'][0, :, :].data

        dataset.close()

        # 执行双线性插值
        precipitation_sta_bilinear = Bilinear_interp(lonSta, latSta, longitude, latitude, precipitation)

        # 将插值结果添加到站点信息 DataFrame
        column_name = f'precipitation_bilinear'
        stations_info[column_name] = precipitation_sta_bilinear

        output_dir = r'E:\satlite\MSWEP\inter2'
        # 构造输出文件名，使用年份和原始nc文件名（不含扩展名）
        output_filename = f"{os.path.splitext(filename)[0]}.xlsx"
        output_path = os.path.join(output_dir, output_filename)

        # 保存到 Excel 文件，每个文件使用年份和原始nc文件名
        try:
            stations_info.to_excel(output_path, sheet_name='Sheet1', engine='openpyxl', index=False)
            print(f"Saved file: {output_path}")
        except PermissionError as e:
            print(f"保存文件时发生权限错误：{e}")
        except Exception as e:
            print(f"保存文件时发生未知错误：{e}")

        # 重置DataFrame以保存下一年的文件
        stations_info.drop(columns=[column_name], inplace=True)

注意：

1.最好先打印一下纬度信息，看一下顺序是从大到小还是从小到大，插值需要纬度是从小到大。如果顺序不是从小到大，则需要翻转纬度

2.有时会出现显示维度不匹配的错误，则需要检查

print("iSta:", iSta, "jSta:", jSta)
print("var shape:", var.shape)

因为可能原始数据不是（lon,lat），可能是（lat,lon）