【python由站点数据插值到网格数据方法对比】

文章【python由站点数据插值到网格数据方法对比】分享给大家，欢迎收藏Python资料网，专注分享技术知识

1、前言

• 气象海洋中空间数据类型有站点数据、格点数据。
• 站点数据空间分布不连续，不利于进行时空分析；有时需要将站点数据插值到网格中。
• 本文基于一个散点数据对比了griddata、IDW、krige、rbf的插值结果。
• interp2d插值，地理数据中的分辨率转换

2、结果对比

2.1 原始散点站位图

• 原始散点站位如下：
  

2.2 griddata插值

from scipy.interpolate import griddata
lat = df_filter['Lat']
lon = df_filter['Lon']
VIS_Min = df_filter['VIS_Min']
olon = np.linspace(lon.min(), lon.max(), 40)
olat = np.linspace(lat.min(), lat.max(), 40)
olon, olat = np.meshgrid(olon, olat)
VIS_Min_grd=griddata((lon,lat),VIS_Min,(olon, olat),method='linear')
plt.figure(figsize=(8,6))
plt.pcolormesh(olon, olat, VIS_Min_grd)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap='viridis',edgecolors='k')

• 可视化代码

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib as mpl

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.colors import LinearSegmentedColormap

import regionmask
import geopandas as gpd
import warnings

warnings.filterwarnings("ignore")
import matplotlib.ticker as ticker
extents = [115.6, 117.2, 39.7, 40.8]
crs = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111, projection=crs)
ax.set_extent(extents, crs)

prov_reader = shpreader.Reader(r'.\bou2_4p.shp')
nineline_reader = shpreader.Reader(r'.\九段线.shp')

ax.add_geometries(prov_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)
ax.add_geometries(nineline_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)

ax.set_xticks(np.arange(extents[0], extents[1] + 0.3, 0.3))
ax.set_yticks(np.arange(extents[2], extents[3] + 0.2, 0.2))
plt.tick_params(axis='both', which='major', labelsize=18)


norm = mcolors.Normalize(vmin=0, vmax=np.round(VIS_Min.max()))
cmap_jet = plt.cm.jet
cmap_gray = plt.cm.gray
gray = cmap_gray(np.linspace(0, 1, 9))
colors = np.vstack((gray[8], cmap_jet(np.linspace(0, 1, 9))))
cmaps = LinearSegmentedColormap.from_list('Combined', colors, N=64)

map = plt.pcolormesh(olon, olat, VIS_Min_grd, cmap=cmaps)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap=cmaps,edgecolors='k',norm=norm)
cb_ax = inset_axes(ax, width="3%", height="100%", loc='lower left', bbox_to_anchor=(1.01, 0., 1, 1),
                       bbox_transform=ax.transAxes, borderpad=0)
cbar = plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmaps), cax=cb_ax, ax=map)
cbar.ax.set_title('能见度', size=18, fontname='SimHei', pad=10)
cbar.ax.tick_params(labelsize=15)
ax.set_xlabel('Lon(°E)', fontsize=18)
ax.set_ylabel('Lat(°N)', fontsize=18)
ax.set_title('原始散点站位', size=18, fontname='SimHei', pad=10)

• 可视化结果
  

2.3 krige插值

• 插值核心代码

from pykrige.ok import OrdinaryKriging
lat = df_filter['Lat']
lon = df_filter['Lon']
VIS_Min = df_filter['VIS_Min']
olon = np.linspace(lon.min(), lon.max(), 40)
olat = np.linspace(lat.min(), lat.max(), 40)

Krin=OrdinaryKriging(lon,lat,VIS_Min,variogram_model='gaussian',
                         coordinates_type="geographic",
                         nlags=1,)
data2,ssl=Krin.execute('grid',olon,olat)

olon, olat = np.meshgrid(olon, olat)
plt.figure(figsize=(8,6))
plt.pcolormesh(olon, olat, data2)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap='viridis',edgecolors='k')

• 插值结果
  

2.4 RBF插值

• 插值核心代码

from scipy.interpolate import Rbf
lat = df_filter['Lat']
lon = df_filter['Lon']
VIS_Min = df_filter['VIS_Min']
olon = np.linspace(lon.min(), lon.max(), 40)
olat = np.linspace(lat.min(), lat.max(), 40)
olon, olat = np.meshgrid(olon, olat)
func1=Rbf(lon,lat,VIS_Min,function='linear')
VIS_Min_rbf=func1(olon,olat)

plt.figure(figsize=(8,6))
plt.pcolormesh(olon, olat, VIS_Min_rbf)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap='viridis',edgecolors='k')

• 插值结果
  

2.5 IDW插值

• 插值核心代码

from math import radians, sin, cos, sqrt, asin
import numpy as np

# degree to km (Haversine method)
def haversine(lon1, lat1, lon2, lat2):
    R = 6372.8  # Earth radius in kilometers
    dLon = radians(lon2 - lon1)
    dLat = radians(lat2 - lat1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
    c = 2*asin(sqrt(a))
    d = R * c
    return d

def IDW(x, y, z, xi, yi):
    lstxyzi = []
    for p in range(len(xi)):
        lstdist = []
        for s in range(len(x)):
            d = (haversine(x[s], y[s], xi[p], yi[p]))
            lstdist.append(d)
        sumsup = list((1 / np.power(lstdist, 2)))
        suminf = np.sum(sumsup)
        sumsup = np.sum(np.array(sumsup) * np.array(z))
        u = sumsup / suminf
        xyzi = [xi[p], yi[p], u]
        lstxyzi.append(xyzi)
    return(lstxyzi)

lat = df_filter['Lat']
lon = df_filter['Lon']
VIS_Min = df_filter['VIS_Min']
olon = np.linspace(lon.min(), lon.max(), 40)
olat = np.linspace(lat.min(), lat.max(), 40)
olon, olat = np.meshgrid(olon, olat)
VIS_Min_idw=IDW(lon,lat,VIS_Min,olon.flatten(), olat.flatten())

VIS_Min_idw= np.array(VIS_Min_idw)[:,2].reshape(olon.shape)
plt.figure(figsize=(8,6))
plt.pcolormesh(olon, olat, VIS_Min_idw)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap='viridis',edgecolors='k')

• 矩阵运算

import numpy as np

def haversine(lon1, lat1, lon2, lat2):
    R = 6372.8  # Earth radius in kilometers
    dLon = np.radians(lon2 - lon1)
    dLat = np.radians(lat2 - lat1)
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    a = np.sin(dLat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dLon / 2) ** 2
    c = 2 * np.arcsin(np.sqrt(a))
    d = R * c
    return d

def IDW(x, y, z, xi, yi):
    xi, yi = np.meshgrid(xi, yi)  # Create grid
    xi = xi.ravel()
    yi = yi.ravel()
    
    dists = haversine(np.tile(x, (len(xi), 1)).T, np.tile(y, (len(yi), 1)).T, xi, yi)
    weights = 1 / dists ** 2
    weights[np.isinf(weights)] = 1e12  # Handle division by zero
    
    weighted_sum = np.sum(weights * z[:, None], axis=0)
    sum_of_weights = np.sum(weights, axis=0)
    
    u = weighted_sum / sum_of_weights
    result = np.column_stack((xi, yi, u))
    
    return result

# 示例数据
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 2, 3, 4])
z = np.array([1, 2, 3, 4, 5])
xi = np.linspace(0, 4, 10)
yi = np.linspace(0, 4, 10)

# 进行IDW插值
interpolated_data = IDW(x, y, z, xi, yi)

print(interpolated_data[:,2])

• 插值结果
  

3、总结

• 插值范围上，除了griddata以外的其它插值方法，都可以插值到站点以外的区域；
• 插值效果上，从本案例的结果上看，我认为反距离权重IDW结果最好，之后为RBF，其次是克里金。
  

4、其它

• 在数据插值过程中，发现了MetPy库下的inverse_distance_to_grid插值函数。用于执行基于反距离权重（Inverse Distance Weighting, IDW）的插值算法，常用于气象学和其他领域的空间数据分析。这个函数允许你将一系列离散观测点的数据插值到一个网格上，从而生成连续的空间场。
• 使用上述数据，调整了相关参数，做了插值，结果并不理想。

from metpy.interpolate import inverse_distance_to_grid
data = inverse_distance_to_grid(lon,lat,VIS_Min, olon, olat, r=0.5)
plt.figure(figsize=(8,6))
plt.pcolormesh(olon, olat, data)
plt.scatter(lon, lat, s=50, c=VIS_Min, cmap='viridis',edgecolors='k')

5、国家站能见度插值结果对比

• 基于全国2000多各国家站逐小时能见度结果进行插值结果对比。
• 对比方法包括griddata、rbf、IDW
• 插值效率：griddata > rbf > IDW
• 插值效果：rbf > griddata > IDW
  

版权声明：本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若内容造成侵权/违法违规/事实不符，请联系邮箱：jacktools123@163.com进行投诉反馈，一经查实，立即删除！