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

本次分享是在上一期的基础上将克里金差值结果进行输出为tif

一、数据介绍

本期使用的数据依然为上一期的所使用的fake数据

二、代码部分

1. 克里金差值部分

克里金差值的核心部分依然是上次所说的OrdinaryKriging函数

from pykrige.ok import OrdinaryKriging

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)

2. tif文件生成部分

tif文件的操作部分我们采用rasterio这个函数来进行，网上很多资料都是根据GDAL这个库来进行的，GDAL这个库操作地理方面的东西确实是鼻祖般的存在，但是对新手不友好，所以今天我们来介绍一下rasterio这个库

import rasterio

如果小伙伴有Anaconda的话直接conda install rasterio就行了，如果没有的话，安装rasterio库需要先安装GDAL这个库。
核心的tif生成函数如下所示：

with rasterio.open(
        'D:\CSDN\克里金插值/测试数据.tif',
        'w',
        driver='GTiff',
        height=len(lat),
        width=len(lon),
        count=1,
        dtype=data.dtype,
        crs='+proj=latlong',
        transform=transform,) as f:
        meta = f.meta
        f.write(data, 1)

三. 分步讲解

1. 库函数引用

import pandas as pd
import rasterio
from rasterio.mask import mask
from rasterio.transform import Affine
from pykrige.ok import OrdinaryKriging
import fiona
import numpy as np

2. 温度数据读取并插值

该部分详细介绍在上一期有讲过

df = pd.read_excel(r'D:\CSDN\克里金插值/假的温度数据.xlsx')
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)

3.transform生成

res = 0.5 #你自己设置的图像分辨率
transform = Affine.translation(lon[0] - res / 2, lat[0] - res / 2) * Affine.scale(res, res) #tif的经纬度生成
# lon 经度，lat 维度 剩下的套用我这句话就可以了

4.tif文件生成

with rasterio.open(
        'D:\CSDN\克里金插值/测试数据.tif',
        'w',
        driver='GTiff', #设置数据类型
        height=len(lat), #维度的格子数量
        width=len(lon), #经度的格子数量
        count=1, #波段数量
        dtype=data.dtype, #数据类型
        crs='+proj=latlong', # 投影
        transform=transform,) as f:
        meta = f.meta
        f.write(data, 1) #想文件中写入数据

用该段代码生成的tif是未经过裁剪的图片

5. tif文件裁剪

# 读取shp文件
with fiona.open(r'D:\CSDN\克里金插值\江苏shp/江苏.shp', "r") as shapefile:
	shapes = [feature["geometry"] for feature in shapefile]
# 载入刚刚输出的测试数据
dsin = rasterio.open(r'D:\CSDN\克里金插值/测试数据.tif')
out_image, _ = mask(dsin, shapes, crop=True, nodata=np.nan) #进行掩膜处理
out_meta = meta #获取tif文件信息
with rasterio.open(r'D:\CSDN\克里金插值/测试数据_裁剪.tif', "w", **out_meta) as dest:
	dest.write(out_image) #输出裁剪好的tif文件

周围黑黑的都是被填充为nan数值了，如果需要做进一步处理的话并不需要在意。