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

文章目录

1、前言
2、结果对比
- 2.1 原始散点站位图
- 2.2 griddata插值
- 2.3 krige插值
- 2.4 RBF插值
- 2.5 IDW插值
3、总结

1、前言

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

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 mplimport cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.colors import LinearSegmentedColormapimport regionmask
import geopandas as gpd
import warningswarnings.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 kilometersdLon = 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)**2c = 2*asin(sqrt(a))d = R * creturn ddef 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 / suminfxyzi = [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')