利用griddata进行插值

全栈程序员站长

发布于 2022-08-24 19:41:02

86400

代码可运行

文章被收录于专栏：全栈程序员必看全栈程序员必看

运行总次数：0

代码可运行

大家好，又见面了，我是你们的朋友全栈君。

利用griddata进行插值

griddata函数讲解
- 第一步：导入相关库
- 第二步：给出插值到的经纬度信息（目标经纬度）
- 第三步：待插值数据
- 第四步：插值
- 汇总成函数
- 结果对比
  - 插值前（10km）
  - 插值后（1km）

因为最近在做算法优化，所以对数据统一性有一定要求，在最近的研究中主要用一个简单的最近邻插值对数据集进行降尺度处理。主要运用到的函数时scipy里面的 griddata

griddata函数讲解

`scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)`

points:2-D ndarray of floats with shape (n, D), 或 length D tuple of 1-D ndarrays with shape (n,).
values:ndarray of float 或 complex, shape (n,)
xi:2-D ndarray of floats with shape (m, D), 或 length D tuple of ndarrays broadcastable to the same shape.
method:{‘linear’, ‘nearest’, ‘cubic’}, 可选参数

nearest:返回最接近插值点的数据点的值 linear:线性插值 cubic:三次样条插值

第一步：导入相关库

import xarray as xr    
from scipy.interpolate import griddata    # 插值函数
import numpy as np

第二步：给出插值到的经纬度信息（目标经纬度）

mask_tmp = xr.open_dataset('G:/China/temperature_max/nc/2000/tmx_2000_1.nc')
# 待插值的标准经纬度
mask_tmp_lon = mask_tmp.lon.values  # (7680,) 一维
mask_tmp_lat = mask_tmp.lat.values  # (4717,) 一维
mask_LON,mask_LAT = np.meshgrid(mask_tmp_lon,mask_tmp_lat)  # (4717, 7680) 生成经纬度网格

第三步：待插值数据

rad = xr.open_dataset('D:/Users/62692/Desktop/rad.nc')
# 辐射数据经纬度
rad_lon = rad.lon.values    # 辐射数据经度 (641,)
rad_lat = rad.lat.values    # 辐射数据纬度 (394,)
rad_LON, rad_LAT = np.meshgrid(rad_lon,rad_lat)   # (394, 641)
rad_LON = rad_LON.ravel().reshape(-1,1)  # 展平 (252554, 1)
rad_LAT = rad_LAT.ravel().reshape(-1,1)  # 展平 (252554, 1)
rad_values = rad['rad'].values     # 需要插值的辐射数据 (394, 641)
points = np.concatenate([rad_LON,rad_LAT],axis = 1)   # (252554, 2)

第四步：插值

# 插值
data = griddata(points, rad_values.ravel(),(mask_LON,mask_LAT),method='nearest')  # 用最近邻插值即可
# rad_values.ravel() (252554,)
# 这里用最邻近主要考虑到辐射数据的连续性变化，对于线性插值或者三次插值并没有多大影响

汇总成函数

''' Created on 1 23, 2022 @author: GongHaixing 将一个文件夹里面所有的nc文件进行插值 '''
def interp2D(maskpath,mask_lon='lon',mask_lat='lat',inputpath='', outputpath='',data_lon='lon',data_lat='lat',variable='',interp_method='nearest',save=True):
    """输入插值目标的相关信息以及需要插值的数据 :maskpath: 需要插值到对应数据的数据路径 :mask_lon: 标准数据的经度名称，比如：x，lon :mask_lat: 标准数据的纬度名称，比如：y，lat :inputpath: 需要做插值处理的nc文件所在的目录 :outputpath: 插值完nc文件保存的路径，注意要是'/' :data_lon: 需要做插值数据经度名称，比如：'x'，'lon' :data_lat: 需要做插值数据经度名称，比如：'y'，'lat' :variable:需要做插值数据变量的名称，比如：'tmp'，'ndvi' :interp_method: griddata的插值方法，比如：'nearest','linear','cubic' :save:是否对文件进行存储 """
    #导入相关库
    import xarray as xr
    import os
    from scipy.interpolate import griddata    # 插值函数
    import numpy as np

    ### 目标插值
    mask_data = xr.open_dataset(maskpath)
    mask_data_lon = mask_data[mask_lon].values
    mask_data_lat = mask_data[mask_lat].values
    mask_LON,mask_LAT = np.meshgrid(mask_data_lon,mask_data_lat)
    mask_LON1 = np.array(mask_LON)
    mask_LAT1 = np.array(mask_LAT)

    ### 插值对象
    os.chdir(inputpath)                            # 给出nc文件所在的目录（路径）
    files = os.listdir()
    savepath = outputpath
    for file in files:
        inputfile = xr.open_dataset(file)
        inputfile_lon = inputfile[data_lon].values                                 # 数据的经纬度
        inputfile_lat = inputfile[data_lat].values
        inputfile_LON, inputfile_LAT = np.meshgrid(inputfile_lon,inputfile_lat)
        inputfile_LON = inputfile_LON.ravel().reshape(-1,1)
        inputfile_LAT = inputfile_LAT.ravel().reshape(-1,1)
        inputfile_values = np.array(inputfile[variable].values,dtype=np.float32)     # 需要插值的数据
        points = np.concatenate([inputfile_LON,inputfile_LAT],axis = 1)
        # 插值
        print('开始对'+file+'进行插值')
        inputfile_interp = griddata(points, inputfile_values.ravel(),(mask_LON1,mask_LAT1),method=interp_method).astype(np.float32)  # 用最近邻插值即可，不然数据会nan
        outfile= xr.Dataset(
                      data_vars = { 
   variable:(('lat','lon'),inputfile_interp)},
                      coords = { 
   
                                'lon':('lon',np.array(mask_data_lon)),
                                'lat':('lat',np.array(mask_data_lat))}
                                )
        if save == True:
            outfile.to_netcdf(outputpath+'/'+file)
            print(file+'已经插值成功，且已经保存到'+outputpath+'路径下')
        else:
            print(file+'已经插值成功,但是我没有保存文件')

from interp2D import *
import xarray as xr
import os
import pandas as pd
from scipy.interpolate import griddata    # 插值函数


maskpath = 'H:/China/temperature_max/nc/2000/tmx_2000_1.nc'
mask_lon='lon'
mask_lat='lat'
inputpath='H:/China/LAI/nc'
outputpath='H:/China/LAI/nc_1km'
data_lon='lon'
data_lat='lat'
variable='LAI'
interp_method='nearest'
save=True

interp2D(maskpath=maskpath,mask_lon='lon',mask_lat='lat',inputpath=inputpath, outputpath=outputpath,data_lon='lon',data_lat='lat',variable=variable,interp_method='nearest',save=True)