点击蓝字，关注我们
Python 格点数据插值到站点：
最邻近插值和双三次插值算法
最邻近插值算法
最邻近插值算法是最简单的插值算法，同时也叫零阶插值法。即选择它所映射位置最近的输入像素的灰度值为结果。对二维图像，是去待采样点周围4个相邻像素点中距离最近的1个点的灰度值作为待采样点的像素值。
下面是将网格数据插值到站点的代码示例：
from pathlib import Pathimport pandas as pdimport numpy as npimport netCDF4 as nc# 读取站点信息stations_info = pd.read_excel(r'D:\ML\grid to stations\stations.xlsx', 'Sheet1', index_col=None, na_values=['NA'])# 读取网格数据dataset = nc.Dataset(r"D:\ML\grid to stations\ERA5.TEST.nc")print(dataset)# 经纬度longitude = dataset.variables['longitude'][:].datalatitude = dataset.variables['latitude'][:].data# 温度t = dataset.variables['t'][:, :, :].data # 获取所有时次的温度数据# 将格点范围内的站点筛选出来lonSta, latSta = stations_info['经度'].to_numpy(), stations_info['纬度'].to_numpy()# 定义获取最临近格点坐标索引的方法def nearest_position(stn_lat, stn_lon, lat2d, lon2d):difflat = stn_lat - lat2ddifflon = stn_lon - lon2drad = np.multiply(difflat, difflat) + np.multiply(difflon, difflon)aa = np.where(rad == np.min(rad))ind = np.squeeze(np.array(aa))return tuple(ind)# 将一维的经纬度数据网格二维化lon2D, lat2D = np.meshgrid(longitude, latitude)# 创建一个 DataFrame 用于存储插值数据t_sta_nearest_df = pd.DataFrame(index=range(len(lonSta)), columns=[f'{hour:02d}h' for hour in range(24)])# 对每个站点进行插值计算for i in range(len(lonSta)):t_nearest = []for t_index in range(t.shape[0]): # 对每个时间点进行操作indexSta = nearest_position(latSta[i], lonSta[i], lat2D, lon2D)jSta, iSta = indexSta[0], indexSta[1]t_nearest.append(t[t_index, jSta, iSta]) # 将当前时间点的结果添加到列表中t_sta_nearest_df.loc[i] = t_nearest# 添加时间标题行t_sta_nearest_df.columns = [f'{hour:02d}h' for hour in range(24)]# 将插值数据添加到站点信息 DataFrame 中stations_info = pd.concat([stations_info, t_sta_nearest_df], axis=1)# 将数据保存为新的xlsx文件stations_info.to_excel('D:/ML/grid to stations/weather_station_data.xlsx', sheet_name='Sheet1', index=False)
双三次插值算法
双三次插值算法(Bicubic interpolation)又称立方卷积插值算法，是对双线性插值的改进，是一种比较复杂的插值方式，它不仅考虑到周围4个像素点灰度值的影像，还考虑到它们灰度值变化率的影像。该算法需要利用待采样附近16个像素点的灰度值作三次插值进行计算。
下面是代码示例
from pathlib import Pathimport pandas as pdimport numpy as npimport netCDF4 as ncfrom scipy.interpolate import griddata# 读取站点信息stations_info = pd.read_excel(r'D:\ML\grid to stations\stations.xlsx', 'Sheet1', index_col=None, na_values=['NA'])# 读取网格数据dataset = nc.Dataset(r"D:\ML\grid to stations\ERA5.TEST.nc")# 经纬度longitude = dataset.variables['longitude'][:].datalatitude = dataset.variables['latitude'][:].data# 温度t = dataset.variables['t'][:, :, :].data # 获取所有时次的温度数据# 将格点范围内的站点筛选出来lonSta, latSta = stations_info['经度'].to_numpy(), stations_info['纬度'].to_numpy()# 将一维的经纬度数据网格化lon2D, lat2D = np.meshgrid(longitude, latitude)# 创建一个 DataFrame 用于存储插值数据t_sta_nearest_df = pd.DataFrame(index=range(len(lonSta)), columns=[f'{hour:02d}h' for hour in range(24)])# 对每个站点进行插值计算for i in range(len(lonSta)):t_nearest = []for t_index in range(t.shape[0]): # 对每个时间点进行操作t_values = t[t_index, :, :].flatten() # 获取当前时间点的所有温度值grid_points = np.vstack((lon2D.flatten(), lat2D.flatten())).T # 网格坐标station_point = np.array([[lonSta[i], latSta[i]]]) # 站点坐标t_interp = griddata(grid_points, t_values, station_point, method='cubic') # 双三次插值t_nearest.append(t_interp[0]) # 将插值结果添加到列表中t_sta_nearest_df.loc[i] = t_nearest# 添加时间标题行t_sta_nearest_df.columns = [f'{hour:02d}h' for hour in range(24)]# 将插值数据添加到站点信息 DataFrame 中stations_info = pd.concat([stations_info, t_sta_nearest_df], axis=1)# 将数据保存为新的xlsx文件stations_info.to_excel('D:/ML/grid to stations/weather_station_data_SSC.xlsx', sheet_name='Sheet1', index=False)
后台私信：第八星系
群内每日更新分享数据
进群请勿回复第八星系以外字词
本文编辑：第八星系欣悦
分享、在看与点赞，至少我要拥有一个吧
