使用griddata进行均匀网格和离散点之间的相互插值

插值操作非常常见，数学思想也很好理解。常见的一维插值很容易实现，相对来说，要实现较快的二维插值，比较难以实现。这里就建议直接使用scipy 的griddata函数。

1 griddata函数介绍

官网介绍

2 离散点插值到均匀网格

def interp2d_station_to_grid(lon,lat,data,loc_range = [18,54,73,135],
                             det_grid = 1 ,method = 'cubic'):
    '''
    func : 将站点数据插值到等经纬度格点
    inputs:
        lon: 站点的经度
        lat: 站点的纬度
        data: 对应经纬度站点的 气象要素值
        loc_range: [lat_min,lat_max,lon_min,lon_max]。站点数据插值到loc_range这个范围
        det_grid: 插值形成的网格空间分辨率
        method: 所选插值方法，默认 0.125
    return:
        
        [lon_grid,lat_grid,data_grid]
    '''
    #step1: 先将 lon,lat,data转换成 n*1 的array数组
    lon = np.array(lon).reshape(-1,1)
    lat = np.array(lat).reshape(-1,1)
    data = np.array(data).reshape(-1,1)
    
    #shape = [n,2]
    points = np.concatenate([lon,lat],axis = 1)
    
    #step2:确定插值区域的经纬度网格
    lat_min = loc_range[0]
    lat_max = loc_range[1]
    lon_min = loc_range[2]
    lon_max = loc_range[3]
    
    lon_grid, lat_grid = np.meshgrid(np.arange(lon_min,lon_max+det_grid,det_grid), 
                       np.arange(lat_min,lat_max+det_grid,det_grid))
    
    #step3:进行网格插值
    grid_data = griddata(points,data,(lon_grid,lat_grid),method = method)
    grid_data = grid_data[:,:,0]
    
    #保证纬度从上到下是递减的
    if lat_grid[0,0]<lat_grid[1,0]:
        lat_grid = lat_grid[-1::-1]
        grid_data = grid_data[-1::-1]
    
    return [lon_grid,lat_grid,grid_data]

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import os
from mpl_toolkits.basemap import Basemap
import numpy.ma as ma

r6_p = get_station_data(f1,file_type = 'r6-p')
scatter_station_on_map(station_lon = r6_p[:,1],station_lat = r6_p[:,2],
                       station_value = r6_p [:,4])

new_data = interp2d_station_to_grid(lon = r6_p[:,1],lat = r6_p[:,2],
                                    data = r6_p[:,4],
#                                    method='linear',
                                    )
contourf_data_on_map(new_data[2],new_data[0],new_data[1])

下面为插值前后的数据类型及其大小

• 原始站点数据。降水量越大，站点颜色越深，小圆圈越大。

• method = ‘linear’

• method = ‘cubic’

可以看到，在点比较少的情况下，不同插值方法，结果相差挺大，但降水中心都预测出来了。

3 均匀网格插值到离散点

在气象上，用得更多的，是将均匀网格的数据插值到观测站点，此时，也可以逆向使用 griddata方法插值；这里就不做图显示了。

def grid_interp_to_station(all_data, station_lon,station_lat ,method = 'cubic'):
    '''
    func: 将等经纬度网格值 插值到 离散站点。使用griddata进行插值
    inputs: 
        all_data,形式为：[grid_lon,grid_lat,data] 即[经度网格，纬度网格，数值网格]
        station_lon: 站点经度
        station_lat: 站点纬度。可以是 单个点，列表或者一维数组
        method: 插值方法,默认使用 cubic
    '''
    station_lon = np.array(station_lon).reshape(-1,1)
    station_lat = np.array(station_lat).reshape(-1,1)
    
    lon = all_data[0].reshape(-1,1)
    lat = all_data[1].reshape(-1,1)
    data = all_data[2].reshape(-1,1)
    
    points = np.concatenate([lon,lat],axis = 1)
    
    station_value = griddata(points,data,(station_lon,station_lat),method=method)
    
    station_value = station_value[:,:,0]
    
    return station_value

4 获取最近邻的Index

def get_nearest_point_index(point_lon_lat,lon_grid,lat_grid):
    '''
    func:获取与给定经纬度值的点最近的等经纬度格点的经纬度index
    inputs:
        point_lon_lat: 给定点的经纬度，eg:[42.353,110.137]
        lon_grid: 经度网格
        lat_grid: 纬度网格
    return:
        index: [index_lat,index_lon]
    '''
    #step1: 获取网格空间分辨率;默认纬度和经度分辨率一致
    det = lon_grid[0,1]-lon_grid[0,0]
    
    #step2: 
    point_lon = point_lon_lat[0]
    point_lat = point_lon_lat[1]
    
    lon_min = np.min(lon_grid)
    lat_min = np.min(lat_grid)
#    lat_max = np.max(lat_grid)
    
    index_lat = round((point_lat-lat_min)/det)
    index_lon = round((point_lon-lon_min)/det)
    
    #由于默认的 lat_max值对应的index为0，因此需要反序
    index_lat = lat_grid.shape[0] -index_lat-1
    
    return [int(index_lat),int(index_lon)]

loc_range = [20,50,100,130]

det_grid = 0.25

lat_min = loc_range[0]
lat_max = loc_range[1]
lon_min = loc_range[2]
lon_max = loc_range[3]

lon_grid, lat_grid = np.meshgrid(np.arange(lon_min,lon_max+det_grid,det_grid), 
                   np.arange(lat_min,lat_max+det_grid,det_grid))

#默认使高纬在第一行
lat_grid = lat_grid[-1::-1]

index = get_nearest_point_index([113.2,30],lon_grid,lat_grid)

print(index)

打印输出的index = [80,53], 我们lon_grid和lat_grid去查找一下，对应的经纬度为[113.25,30] ,

刚好位置对上！done!