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

使用griddata进行均匀网格和离散点之间的相互插值文章目录1griddata函数介绍2离散点插值到均匀网格3均匀网格插值到离散点4获取最近邻的Index插值操作非常常见,数学思想也很好理解。常见的一维插值很容易实现,相对来说,要实现较快的二维插值,比较难以实现。这里就建议直接使用scipy的griddata函数。1griddata函数介绍官网介绍2离散点插值到均匀网格definterp2d_station_to_gri…

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

插值操作非常常见,数学思想也很好理解。常见的一维插值很容易实现,相对来说,要实现较快的二维插值,比较难以实现。这里就建议直接使用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!

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/143986.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号