Почему ошибка ValueError?
Прежде всего, scipy.interpolate.interpn требует, чтобы interp_points.shape[-1] было таким же, как число измерений в вашей задаче. Вот почему вы получаете ValueError из своего фрагмента кода - ваш interp_points имеет 92476 как n_dims, что конфликтует с фактическим количеством димсов (4).
Быстрое исправление
Вы можете исправить этот фрагмент кода, просто изменив порядок своих операций. Вы пытаетесь сжать слишком рано - если вы сжимаете после интерп:
points = (solzen, aod, index, wave) # 7, 5, 92476, 140
mg = np.meshgrid(solzen0, aod0, index0, wave0) # 4, 1, 1, 92476, 70
interp_points = np.moveaxis(mg, 0, -1) # 1, 1, 92476, 70, 4
result_presqueeze = interpn(points,
skyrad0, interp_points) # 1, 1, 92476, 70
result = np.squeeze(result_presqueeze,
axis=(0,1)) # 92476, 70
Я заменил interp_mesh на mg здесь и удалил np.array (это не требуется, поскольку np.meshgrid возвращает объект ndarray).
Комментарии по производительности
Я думаю, что ваш фрагмент кода в порядке, однако вы можете sh для использования xarray, если вы обрабатываете помеченные данные, как это:
- чуть более читабельно, чем немаркированные
numpy массивы - также может обрабатывать некоторые из фоновая работа с использованием dask (полезно, если вы просматриваете большие объемы данных в 6D)
Обновление : Ой! Это должно было быть .interp, а не .loc. Фрагмент кода ниже работал, потому что точки данных были фактически исходными точками данных. В качестве предупреждения для других:
from scipy.interpolate import interpn
import numpy as np
from xarray import DataArray
# Define the data space in the 4D skyrad0 array
solzen = np.arange(0,70,10) # 7
aod = np.arange(0,0.25,0.05) # 5
index = np.arange(1,92477,1) # 92476
wave = np.arange(350,1050,5) # 140
# Simulated skyrad for the values above
skyrad0 = np.random.rand(
solzen.size,aod.size,index.size,wave.size) # 7, 5, 92476, 140
# Data space for desired output values of skyrad
# with interpolation between input data space
solzen0 = 30 # 1
aod0 = 0.1 # 1
index0 = index # 92476
wave0 = np.arange(350,1050,10) # 70
def slow():
points = (solzen, aod, index, wave) # 7, 5, 92476, 140
mg = np.meshgrid(solzen0, aod0, index0, wave0) # 4, 1, 1, 92476, 70
interp_points = np.moveaxis(mg, 0, -1) # 1, 1, 92476, 70, 4
result_presqueeze = interpn(points,
skyrad0, interp_points) # 1, 1, 92476, 70
result = np.squeeze(result_presqueeze,
axis=(0,1)) # 92476, 70
return result
# This function uses .loc instead of .interp!
"""
def fast():
da = DataArray(name='skyrad0',
data=skyrad0,
dims=['solzen','aod','index','wave'],
coords=[solzen, aod, index, wave])
result = da.loc[solzen0, aod0, index0, wave0].squeeze()
return result
"""
Сделав несколько изменений в обновленном фрагменте кода, заданном OP:
import numpy as np
import xarray as xr
from scipy.interpolate import interpn
azimuth = np.arange(0, 185, 5) # 37
senzen = np.arange(0, 185, 5) # 37
#wave = np.arange(350,1050,5) # 140
wave = np.asarray([350, 360, 370, 380, 390, 410, 440, 470, 510,
550, 610, 670, 750, 865, 1040, 1240, 1640, 2250]) # 18
solzen = np.arange(0,65,5) # 13
aod = np.arange(0,0.55,0.05) # 11
wind = np.arange(0, 20, 5) # 4
coords = [azimuth, senzen, wave, solzen, aod, wind]
azimuth0 = 135 # 1
senzen0 = 140 # 1
wave0 = np.arange(350,1010,10) # 66
solzen0 = 30 # 1
aod0 = 0.1 # 1
wind0 = 10 # 1
interp_coords = [azimuth0, senzen0, wave0, solzen0, aod0, wind0]
# Simulated rad_boa
rad_boa = np.random.rand(
*map(lambda x: x.size, coords)) # 37, 37, 140/18, 13, 11, 4
def slow():
mg = np.meshgrid(*interp_coords)
interp_points = np.moveaxis(mg, 0, -1)
result_presqueeze = interpn(coords,
rad_boa, interp_points)
result = np.squeeze(result_presqueeze)
return result
def fast():
da = xr.DataArray(name='Radiance_BOA',
data=rad_boa,
dims=['azimuth','senzen','wave','solzen','aod','wind'],
coords=coords)
interp_dict = dict(zip(da.dims, interp_coords))
rad_inc_scaXR = da.interp(**interp_dict).squeeze()
return rad_inc_scaXR
Это довольно быстро:
>>> %timeit slow()
2.09 ms ± 85.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
>>> %timeit fast()
343 ms ± 6.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>>> np.array_equal(slow(),fast())
True
Подробнее о xarray интерполяции можно узнать здесь . Экземпляры набора данных имеют очень похожий синтаксис.
Также возможно изменить метод интерполяции по желанию (возможно, можно sh предоставить ключевой аргумент method='nearest' - .interp для дискретной задачи интерполяции) .
Более продвинутые вещи
Если вы хотите реализовать что-то более продвинутое, я бы порекомендовал, возможно, использовать одну из реализаций MARS (многомерных адаптивных сплайнов регрессии) , Он находится где-то между стандартной регрессией и интерполяцией и работает для многомерных случаев. В Python 3 ваша лучшая ставка - pyearth.