Python站点数据插值与网格数据方法对比

文章目录

1、前言

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')

插值结果

3、总结

插值效果上，从本案例的结果上看，我认为反距离权重IDW结果最好，之后为RBF，其次是克里金。

4、其它

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、国家站能见度插值结果对比

插值效果：rbf > griddata > IDW

作者：木叶清风666