python中带有网格数据的线性插值

Question

我有一个网格天气数据集，其维度为 33 X 77 X 77。第一个维度是时间，休息时间分别是纬度和经度。我需要每次将数据插值（线性或最近邻）到不同的点（纬度和经度）并将其写入 csv 文件。我已经使用了 scipy 的 interp2d 函数，并且它在一个时间步骤中是成功的。由于我有很多地方，我不想随着时间的推移而循环。

下面显示的是我编写的一段代码，有人可以提出更好的方法来完成任务吗？

import sys ; import numpy as np ; import scipy as sp ; from scipy.interpolate import interp2d ;import datetime ; import time ; import pygrib as pg ; 
grb_f=pg.open('20150331/gfs.20150331.grb2')  lat=tmp[0].data(lat1=4,lat2=42,lon1=64,lon2=102)[1] ; lat=lat[:,0]; 
lon=tmp[0].data(lat1=4,lat2=42,lon1=64,lon2=102)[2] ; lon=lon[0,:] ; 
temp=np.empty((0,lon.shape[0]))
for i in range(0,tmp.shape[0]):
    dat=tmp[i].data(lat1=4,lat2=42,lon1=64,lon2=102)
    temp=np.concatenate([temp,dat[0]-273.15],axis=0)
temp1=temp.reshape(tmp.shape[0],lat.shape[0],lon.shape[0])
x=77 ; y=28 #(many points) 
f=interp2d(lon,lat, temp1[0,:,:],kind='linear',copy=False,bounds_error=True ) ; Z=f(x,y)  

编辑::

我没有制作 3D 矩阵，而是将数据垂直添加并制作了大小为 2541 X 77 的数据矩阵以及大小为 2541 X 1 的 lat 和 lon。interp2d 函数给出了无效长度错误。

f=interp2d(lon,lat, temp1[0,:,:],kind='linear',copy=False,bounds_error=True )

“非矩形网格的输入 z 长度无效”)

ValueError：非矩形网格的输入 z 长度无效

我的 x,y,z 矩阵的长度相同 (2541,2541,2541)。那为什么它会抛出一个错误？ 谁能解释一下？我们将非常感谢您的帮助。

Answer 1

如果每次的纬度和经度都相同，您是否可以使用切片和手动插值来实现。所以如果你想要一个 lat = 4.875, lon = 8.4 的一维数组（显然你需要缩放以匹配你的实际间距）

b = a[:,4:6, 8:10]
c = ((b[:,0,0] * 0.125 + b[:,0,1] * 0.875) * 0.6 + ((b[:,1,0] * 0.125 + b[:,1,1] * 0.875) * 0.4)

显然你可以在一行中完成所有操作，但它会更丑

编辑以允许在每个时间段使用变量 lat 和 lon。

lat = np.linspace(55.0, 75.0, 33)
lon = np.linspace(1.0, 25.0, 33)
data = np.linspace(18.0, 25.0, 33 * 77 * 77).reshape(33, 77, 77)

# NB for simplicity I map 0-360 and 0-180 rather than -180+180
# also need to ensure values on grid lines or edges work ok
lat_frac = lat * 77.0 / 360.0
lat_fr = np.floor(lat_frac).astype(int)
lat_to = lat_fr + 1
lat_frac -= lat_fr

lon_frac = lon * 77.0 / 180.0
lon_fr = np.floor(lon_frac).astype(int)
lon_to = lon_fr + 1
lon_frac -= lon_fr

data_interp = ((data[:,lat_fr,lon_fr] * (1.0 - lat_frac) +
                data[:,lat_fr,lon_to] * lat_frac) * (1.0 - lon_frac) +
               (data[:,lat_to,lon_fr] * (1.0 - lat_frac) +
                data[:,lat_to,lon_to] * lat_frac) * lon_frac)

Answer 2

如果您想创建一次“插值器”对象，并使用它来顺序查询您需要的特定点，您可以在scipy.interpolate.Rbf 模块中抢夺：

“n维分散数据的径向基函数逼近/插值类。”

如果您调整时间维度和空间维度之间的比率，n-dimensional 将适用于您的数据，scattered 意味着您也可以将其用于常规/统一数据。

Answer 3

使用RedBlackPy 处理时间序列非常容易。

 import datetime as dt
 import redblackpy as rb

 index = [dt.date(2018,1,1), dt.date(2018,1,3), dt.date(2018,1,5)]
 lat = [10.0, 30.0, 50.0]
 # create Series object
 lat_series = rb.Series(index=index, values=lat, dtype='float32',
                        interpolate='linear')
 # Now you can access at any key using linear interpolation
 # Interpolation does not create new items in Series
 # It uses neighbours to calculate value inplace when you call getitem
 print(lat_series[dt.date(2018,1,2)]) #prints 20

因此，如果您只想将插值写入 csv 文件，您可以遍历所需键的列表并调用系列对象的 getitem，然后将值放入文件：

 # generator for dates range
 def date_range(start, stop, step=dt.timedelta(1)):

     it = start - step

     while it < step:
         it += step
         yield it

 #------------------------------------------------
 # create list for keeping output strings
 out_data = []    
 # create output file
 out_file = open('data.csv', 'w')
 # add head for output table
 out_data.append('Time,Lat\n')

 for date in date_range(dt.date(2018,1,1), dt.date(2018,1,5)):
     out_data.append( '{:},{:}\n'.format(date, lat_series[date]) )

 # write output Series
 out_file.writelines(out_data)
 out_file.close()

您可以通过同样的方式将 Lon 数据添加到您的处理中。