数据处理与可视化 | 站点插值格点+空间区域掩膜

郭好奇同学

发布于 2021-02-12 18:09:00

2.1K0

发布于 2021-02-12 18:09:00

文章被收录于专栏：好奇心Log

前两天写了插值+空间掩膜的推文，不过因为数据问题删除了。后台很多朋友留言说有需要，还是想学习一下，因此自己造了个数据再把这篇文章推一遍。

站点->格点

空间数据类型有站点数据、格点数据。站点数据位置分布不规则，不利于可视化与分析，所以有时需要把站点数据转化成规则的、连续的格点数据。常用的插值方法有克里金插值、径向基插值、反向权重插值......反正还挺多的，今天打算介绍一下克里金插值。

克里金插值

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

像这种常用的算法一般都是有前人做过的，在github上一搜就能查到。对于克里金插值可以直接调用pykrige包进行Kriging插值计算。

代码

# coding: utf-8
import math
import numpy as np
import pandas as pd
import geopandas as gpd
import xarray as xr
import salem
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from pykrige.ok import OrdinaryKriging
import os

def create_map(fig):
    shp_path = r'/China_shp/'
    # --创建画图空间
    proj = ccrs.LambertConformal(central_longitude=110, central_latitude=35, 
                             standard_parallels=(30, 60))  # 创建坐标系

    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
    # --设置shp
    provinces = cfeat.ShapelyFeature(
        Reader(shp_path + 'province.shp').geometries(),
        ccrs.PlateCarree(), 
        edgecolor='k',
        facecolor='none',
    )
    ax.add_feature(provinces, linewidth=0.6, zorder=2)
    # --设置范围
    ax.set_extent([80, 135, 15, 55], crs=ccrs.PlateCarree())
    # --关闭边框
    ax.spines['geo'].set_visible(False)
    # --设置刻度
    ax.tick_params(
        axis='both', 
        labelsize=5, 
        direction='out',
        labelbottom=False, 
        labeltop=False, 
        labelleft=False, 
        labelright=False
    )
    # --设置小地图
    left, bottom, width, height = 0.675, 0.13, 0.23, 0.27
    ax2 = fig.add_axes(
        [left, bottom, width, height], 
        projection=ccrs.PlateCarree()
    )
    provinces = cfeat.ShapelyFeature(
        Reader(shp_path + 'province.shp').geometries(),
        ccrs.PlateCarree(), 
        edgecolor='k',
        facecolor='gray',
    )
    ax2.add_feature(provinces, linewidth=0.6, zorder=2)
    ax2.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())    
    return ax


def krige_interp(df):
    #插值成0.25*0.25
    grid_lon = np.arange(float(math.floor(df.lon.min())), float(math.ceil(df.lon.max())), 0.25)
    grid_lat = np.arange(float(math.floor(df.lat.min())), float(math.ceil(df.lat.max())), 0.25)
    lons = df["lon"].values
    lats = df["lat"].values
    data = df["data"].values
    
    #调用克里金插值函数
    OK = OrdinaryKriging(lons, lats, data,
                         variogram_model='gaussian',
                         coordinates_type="geographic",
                         nlags=6,)
    z, ss = OK.execute('grid', grid_lon, grid_lat)
    #转换成网格
    xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
    #将插值网格数据整理
    df_grid = pd.DataFrame(dict(lon=xgrid.flatten(), lat=ygrid.flatten()))
    #添加插值结果
    df_grid['data'] = z.flatten()
    df_grid = df_grid.set_index(['lat', 'lon'])
    ds = xr.Dataset.from_dataframe(df_grid).data
    return ds


def mask_region(ds, region='china'):
    shp_path = r'/China_shp'
    shp = gpd.read_file(shp_path + 'province.shp').set_index('省')
    if region != 'china':
        shp = shp.loc[region, :]
    ds_mask = ds.salem.roi(shape=shp)
    return ds_mask


def draw_fig(ds, levels, cmap, outpath='.'):
    fig = plt.figure(figsize=(10, 7), dpi=130)  # 创建页面
    
    im = ds.plot.contourf(
        ax=create_map(fig), 
        cmap=cmap, 
        transform=ccrs.PlateCarree(),
        levels=levels,
        add_colorbar=False
    )
    cbar_ax = fig.add_axes([0.13, 0.17, 0.3, 0.012])
    cb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
    plt.savefig(f'{outpath}{os.sep}kriging_mask_chinamap.png')
    
if __name__ == "__main__":
    data_file = r'/testData.csv'
    levels = np.arange(116, 123, 0.5)
    cmap = 'Spectral_r'

    df = pd.read_csv(data_file, header=0, names=['lon', 'lat', 'data'])
    ds = krige_interp(df)
    ds_mask = mask_region(ds)
    draw_fig(ds_mask, levels, cmap,)

站点的密集程度对插值效果影响非常之大，由于中西部数据量比较少，所以插值结果也是没眼看。

也可以进一步调整，只留下部分省份，其余地区掩盖掉。

if __name__ == "__main__":
    data_file = r'./data.csv'
    levels = np.arange(116, 123, 0.5)
    cmap = 'Spectral_r'

    df = pd.read_csv(data_file, header=0, names=['lon', 'lat', 'data'])
    ds = krige_interp(df)
    ds_mask = mask_region(ds, region=['北京市', '天津市', '河北省', '山西省', '山东省', '辽宁省'])
    draw_fig(ds_mask, levels, cmap,)