Python气象数据处理与绘图：克里金(Kriging)插值与可视化

气海无涯 2021-11-25

3332

克里金法（Kriging）是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计（Best Linear Unbiased Prediction, BLUP），因此在地统计学中也被称为空间最优无偏估计器(spatial BLUP)。

1、安装模块

1!pip install PyKrige
2!pip install plotnine
3!pip install openpyxl

复制

2、导入模块

 1import pandas as pd
 2from pykrige.ok import OrdinaryKriging
 3import plotnine
 4from plotnine import *
 5import geopandas as gpd
 6import shapefile
 7import matplotlib.pyplot as plt
 8import cartopy.crs as ccrs
 9import cartopy.io.shapereader as shpreader
10import cmaps
11from matplotlib.path import Path
12from matplotlib.patches import PathPatch

复制

3、读取数据

1df = pd.read_excel('/home/mw/input/meiyu6520/meiyu_sh_2020.xlsx')
2df

复制

1# 读取站点经度
2lons = df['lon']
3# 读取站点纬度
4lats = df['lat']
5# 读取梅雨量数据
6data = df['meiyu']
7# 生成经纬度网格点
8grid_lon = np.linspace(120.8, 122.1,1300)
9grid_lat = np.linspace(30.6, 31.9,1300)

复制

4、克里金(Kriging)插值

1OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian',nlags=6)
2z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
3z1.shape

复制

输出：

1(1300, 1300)

复制

转换成网格

1xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)

复制

将插值网格数据整理

1df_grid = pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))

复制

添加插值结果

1df_grid["Krig_gaussian"] = z1.flatten()
2df_grid

复制

输出：

 1long    lat Krig_gaussian
 20    121.000000  30.5    541.076134
 31    121.020408  30.5    540.997891
 42    121.040816  30.5    540.770766
 53    121.061224  30.5    540.381399
 64    121.081633  30.5    539.818047
 7...    ... ... ...
 82495    121.918367  31.5    567.391506
 92496    121.938776  31.5    562.730899
102497    121.959184  31.5    558.418169
112498    121.979592  31.5    554.460195
122499    122.000000  31.5    550.857796
132500 rows × 3 columns

复制

5、读取上海的行政区划

1sh = gpd.read_file('/home/mw/input/meiyu6520/Shanghai.shp')
2sh

复制

输出：

 1City    District    Province    Code    geometry
 20    上海市 普陀区 上海市 310107  POLYGON ((121.35622 31.23362, 121.35418 31.237...
 31    上海市 宝山区 上海市 230506  POLYGON ((121.48552 31.31156, 121.48541 31.311...
 42    上海市 崇明区 上海市 310151  MULTIPOLYGON (((121.87022 31.29554, 121.86596 ...
 53    上海市 奉贤区 上海市 310120  POLYGON ((121.56443 30.99643, 121.57047 30.998...
 64    上海市 虹口区 上海市 310109  POLYGON ((121.46828 31.31520, 121.46831 31.316...
 75    上海市 黄浦区 上海市 310101  POLYGON ((121.48781 31.24419, 121.49483 31.242...
 86    上海市 嘉定区 上海市 310114  POLYGON ((121.33689 31.29506, 121.33650 31.294...
 97    上海市 金山区 上海市 310116  POLYGON ((121.00206 30.95104, 121.00764 30.947...
108    上海市 静安区 上海市 310106  POLYGON ((121.46808 31.32032, 121.46809 31.320...
119    上海市 闵行区 上海市 310112  POLYGON ((121.33689 31.23674, 121.33835 31.237...
1210    上海市 浦东新区    上海市 310115  MULTIPOLYGON (((121.97077 31.15756, 121.96568 ...
1311    上海市 青浦区 上海市 310118  POLYGON ((121.31900 31.15867, 121.31953 31.157...
1412    上海市 松江区 上海市 310117  POLYGON ((121.02906 30.94388, 121.02804 30.943...
1513    上海市 徐汇区 上海市 310104  POLYGON ((121.45800 31.21929, 121.45807 31.219...
1614    上海市 杨浦区 上海市 310110  POLYGON ((121.52288 31.34289, 121.52549 31.346...
1715    上海市 长宁区 上海市 310105  POLYGON ((121.43938 31.21444, 121.43946 31.214...

复制

绘图：

1sh.plot()

复制

6、插值结果可视化

 1fig = plt.figure(figsize=(12, 8))
 2ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
 3ax.set_extent([120.8, 122.1, 30.6, 31.9], crs=ccrs.PlateCarree())
 4
 5province = shpreader.Reader('/home/mw/input/meiyu6520/Shanghai.shp')
 6ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')
 7
 8cf = ax.contourf(xgrid, ygrid, z1, levels=np.linspace(0,800,21), cmap=cmaps.MPL_rainbow, transform=ccrs.PlateCarree())
 9
10def shp2clip(originfig, ax, shpfile):
11    sf = shapefile.Reader(shpfile)
12    vertices = []
13    codes = []
14    for shape_rec in sf.shapeRecords():
15        pts = shape_rec.shape.points
16        prt = list(shape_rec.shape.parts) + [len(pts)]
17        for i in range(len(prt) - 1):
18            for j in range(prt[i], prt[i + 1]):
19                vertices.append((pts[j][0], pts[j][1]))
20            codes += [Path.MOVETO]
21            codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
22            codes += [Path.CLOSEPOLY]
23        clip = Path(vertices, codes)
24        clip = PathPatch(clip, transform=ax.transData)
25    for contour in originfig.collections:
26        contour.set_clip_path(clip)
27    return contour
28
29shp2clip(cf, ax, '/home/mw/input/meiyu6520/Shanghai.shp')
30
31cb = plt.colorbar(cf)
32cb.set_label('Rainfall(mm)',fontsize=15)
33
34plt.show()
复制