用Python实现插值，并根据地图裁剪进行区域裁剪

一、数据介绍

1.台站介绍

本次教程以江苏省为例，将简单介绍一下如何利用台站数据进行插值为面数据，台站信息百度上有很多，我也是随便找了一个下载的。

2.假数据生成

因为是教程，所以我在这里利用Excel的randbetween函数随机生成了一列数据，如下所示：

二、代码部分

1.库函数准备

import pandas as pd
from pykrige.ok import OrdinaryKriging
import shapefile,cmaps
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

这里特别说明一下shapefile这个库在安装的时候需要使用的命令是 pip install pyshp千万不要用错哦，而且这个库是基于GDAL的，如果大家有Anaconda的话还是推荐大家直接使用conda

2. 核心函数OrdinaryKriging介绍

OrdinaryKriging 函数就是我们插值的核心函数，这个函数主要输入变量为：

OrdinaryKriging(经度, 纬度, 待插值变量, 变异函数类型(我在这里用的是：gaussian))

更为详细的参数信息可以点击这里，去官网荡一下！！！

3.数据读取

读取Excel文件自然还是使用国宝(pandas)来读取了，我们所需要的信息也就只有：经度、纬度、温度三列即可。

df = pd.read_excel(r'D:\CSDN\克里金插值/假的温度数据.xlsx')
print(df)

接下来就是插值环节，插值第一件事就是需要知道插值的地理范围，江苏省地理范围是：东经116至122，北纬30至36。

Krin = OrdinaryKriging(df['经度'], df['纬度'], df['温度'], variogram_model='gaussian')
lon = np.arange(116, 122.1, 0.5) # 经度,分辨率为0.5
lat = np.arange(30, 36.1, 0.5) # 纬度，分辨率为0.5
data, ssl = Krin.execute('grid', lon, lat)
# print(Krin.execute('grid', grid_lon, grid_lat))
# print(data)

是的这就是插值环节，很简单的，如果我我们想把插值后的结果再输出为Excel这样的一个经纬度对应一个温度值的话，那么我们就需要吧经纬度和data进行匹配，代码如下：

x_grid, y_grid = np.meshgrid(lon,lat)
df_grid =pd.DataFrame(dict(long=x_grid.flatten(),lat=y_grid .flatten()))
df_grid["Krig_gaussian"] = z1.flatten()

插值的结果看起来也还不错…

4.出图

当然，我作为一名头秃的气象人肯定是要出一张等值线填充图了，出图么，就介绍一下，如何得到江苏省（裁剪）的温度分布图。

plt.contourf(lon,lat,data)
plt.show()

这是最基本的出图，但是！！！！！这种图鬼知道是哪，而且一点也不美观，接下来我们就加上亿点点修饰。

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([116, 122.1, 30, 35.6], crs=ccrs.PlateCarree())

province = shpreader.Reader(r'D:\CSDN\克里金插值\江苏shp/江苏.shp')
c = ax.contourf(lon,lat,data,cmap='coolwarm')
ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5, edgecolor='k',
                facecolor='none')
plt.colorbar(c)

结果了简单的美化，现在有一点点ArcGis出图的那个味道了，接下来一步就是，按照shp文件把地图裁剪出来。

利用shp裁剪地图的核心函数，该函数输入项有你的contourf函数输出，坐标轴ax，以及shp文件路径。

def shp2clip(originfig, ax, shpfile):
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i + 1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)
    for contour in originfig.collections:
        contour.set_clip_path(clip)
    return contour

此处完整代码如下：

plt.rcParams['savefig.dpi'] = 300 # 图片像素
plt.rcParams['figure.dpi'] = 300 # 分辨率
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([116, 122.1, 30, 35.6], crs=ccrs.PlateCarree())
province = shpreader.Reader(r'D:\CSDN\克里金插值\江苏shp/江苏.shp')
c = ax.contourf(lon,lat,data,cmap='coolwarm')
ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5, edgecolor='k',
                facecolor='none')

from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter  # 刻度格式
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))  
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(116,122.1,1),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(30,35.7,1),crs=ccrs.PlateCarree())

shp2clip(c, ax, r'D:\CSDN\克里金插值\江苏shp/江苏.shp')
cb = plt.colorbar(c)
cb.set_label('克里金温度插值', fontsize=15)

plt.show()

终于有GIS出图那个味儿了

三、完整代码奉上

import pandas as pd
from pykrige.ok import OrdinaryKriging
import shapefile,cmaps
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

plt.rcParams['savefig.dpi'] = 300 # 图片像素
plt.rcParams['figure.dpi'] = 300 # 分辨率
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

def shp2clip(originfig, ax, shpfile):
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i + 1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)
    for contour in originfig.collections:
        contour.set_clip_path(clip)
    return contour

df = pd.read_excel(r'D:\CSDN\克里金插值/假的温度数据.xlsx')
#print(df)
Krin = OrdinaryKriging(df['经度'], df['纬度'], df['温度'], variogram_model='gaussian')
lon = np.arange(116, 123.1, 0.5) # 经度,分辨率为0.5
lat = np.arange(30, 36.1, 0.5) # 纬度，分辨率为0.5
data, ssl = Krin.execute('grid', lon, lat)
# print(Krin.execute('grid', grid_lon, grid_lat))
# print(data)

x_grid, y_grid = np.meshgrid(lon,lat)
df_grid =pd.DataFrame(dict(long=x_grid.flatten(),lat=y_grid .flatten()))
df_grid["温度"] = data.flatten()
#print(df_grid)

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([116, 122.1, 30, 35.6], crs=ccrs.PlateCarree())

province = shpreader.Reader(r'D:\CSDN\克里金插值\江苏shp/江苏.shp')
c = ax.contourf(lon,lat,data,cmap='coolwarm')
ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5, edgecolor='k',
                facecolor='none')

ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))  
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(116,122.1,1),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(30,35.7,1),crs=ccrs.PlateCarree())

shp2clip(c, ax, r'D:\CSDN\克里金插值\江苏shp/江苏.shp')
cb = plt.colorbar(c)
cb.set_label('克里金温度插值', fontsize=15)
plt.show()

四、下期预告

划重点！！

下期介绍如何将点插值为面，并且保存为tif文件(含坐标系)