大家好,又见面了,我是你们的朋友全栈君。
插值操作非常常见,数学思想也很好理解。常见的一维插值很容易实现,相对来说,要实现较快的二维插值,比较难以实现。这里就建议直接使用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账号...