Python分析之3 种空间插值方法

插值是一个非常常见的数学概念，不仅数据科学家使用它，而且各个领域的人们也使用它。然而，在处理地理空间数据时，插值变得更加复杂，因为您需要基于几个通常稀疏的观测值创建代表性网格。

在深入研究地理空间部分之前，让我们简要回顾一下线性插值。

为了演示的目的，我将使用正则多项式函数：

def F(x):
  return -2*x**3+x**2+2.1

x = np.arange(-5,5, 0.1)
y = F(x)

现在我们可以随机采样几个点 [-4.2, 0, 2.5] 并将它们连接在一起：

这称为线性插值，因为函数在每个区间都用一条直线来近似，现在，只知道函数在 3 个点的值，我们就可以找到区间 [-4.2;2.5] 内的值。

还有许多其他方法，它们具有更高的精度，但它们背后的想法是相同的：找到至少两个已知点之间的函数值。

现在是时候进入地理空间部分了。在本教程中，我们的目标是对NOAA提供的瑞士各地气象站点测量的每日平均气温进行空间插值。预期结果是温度网格，单元格分辨率为 0.1°。

首先，我们需要获取瑞士的行政边界并使用 geopandas 将其可视化：

import geopandas as gdp

shape = gpd.read_file('gadm41_CHE_0.shp')
shape.plot()

 

现在让我们绘制温度观测值并将其与国家形状叠加。为此，我们将气象数据加载到常规 pandas 数据框中，然后将其转换为 geopandas 数据框，并将坐标转换为形状点：

import pandas as pd
from shapely.geometry import Point

df = pd.read_csv('3639866.csv')

points = list()
for i in range(len(df)):
  point = Point(df.loc[i, 'LONGITUDE'], df.loc[i, 'LATITUDE'])
  points.append(point)

gdf = gpd.GeoDataFrame(geometry=points).set_crs(shape.crs)

完成此操作后，我们可以使用matplotlib轻松覆盖两个数据帧。

fig, ax = plt.subplots(figsize=(16,9))
shape.plot(ax=ax, color='black')
gdf.plot(ax=ax, color='r', markersize=85)
plt.show()

为了可视化我们的任务，让我们创建用于插值的网格并将其与上面的地图重叠：

import matplotlib.ticker as mticker
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

LAT, LON = np.arange(45.75, 48, 0.1), np.arange(6, 10.81, 0.1)

fig, ax = plt