Сохранение данных перед передачей их в процедуру интерполяции

#python #numpy #scipy #interpolation

#python #numpy #scipy #интерполяция

Вопрос:

Я начну с извинений, потому что я совершенно новичок в Python (здесь парень из Fortran) и учился на лету. В результате, вероятно, в моих знаниях есть несколько довольно вопиющих пробелов, которые могут быть очевидны после прочтения моей текущей дилеммы.

У меня есть некоторые данные, которые необходимо будет записать в файл, где их затем можно будет прочитать с помощью алгоритма интерполяции. На данный момент этот алгоритм интерполяции, вероятно, будет RectBivariateSpline из SciPi. Хотя интервалы X и Y не совпадают, они являются регулярными, так что это казалось бы идеальным.

Данные обычно имеют вид

X1,Y1,F(X1,Y1)

X1,Y2,F(X1,Y2)

X1,Y3,F(X1,Y3)

X2,Y1,F(X2,Y1)

X2,Y2,F(X2,Y2)

X2,Y3,F(X2,Y3)

и т.д…

В этом случае F(X, Y) не является явной математической функцией, а скорее точкой данных для физической величины в X, Y.

Данные считываются из источника данных, который я в принципе не могу контролировать, используя

  data_array = np.loadtxt(Path/DataFile, dtype = Float, delimiter = ";", usecols = #) 
  

В этом случае возможно несколько разных столбцов с разными данными, но все они зависят от X и Y. Существует отдельный файл, содержащий информацию о диапазоне и размере шага величин X и Y, которые я считываю и сохраняю в массиве. По крайней мере, я вполне уверен, что это массив, а не список, поскольку я где-то читал, что np.loadtxt и np.genfromtxt генерируют массивы numpy, а не списки ванильного python.

Я перебрал множество идей о том, как оптимально хранить эти данные для передачи на другой компьютер с помощью простой процедуры интерполяции, но мне не помешал бы один совет. Мне впервые пришло в голову использовать

 ArrayExample = np.Empty(xRange,yRange) 
For n in (xRange) 
    For m in (yRange) 
        ArrayExample[n,m] = F(X,Y)
  

Однако это ничего не делает для сохранения фактического количества X или Y, связанного со значением в массиве, а они необходимы для интерполяции, и это определенно необходимо для построения графика.

И тогда мне пришло в голову, что, поскольку у меня есть значения X и Y в форме, которую я могу довольно легко прочитать, я мог бы сделать что-то следующим образом. Где xvalues и yvalues — это массивы, содержащие фактические значения X и Y.

 ArrayExample = np.Empty(xRange,yRange,1) 
For n in (xRange) 
    For m in (yRange) 
        ArrayExample[n,m,1] = (xvalues[n],yvalues[m],F(X,Y))
  

После чего мне приходит в голову сохранить ArrayExample как файл pickle, который будет перемещен куда угодно, а затем введен с помощью pickle на другом конце.

Однако, как только у меня это есть, я действительно не знаю, как заставить RectBivariateSpline принимать данные. Я пытался прочитать документацию на сайте scipy и поискать в Google, но все, что я нашел до сих пор, было довольно бесполезным. Если у кого-нибудь есть хорошие примеры того, как его использовать, это было бы очень полезно.

Мы будем признательны за любые советы, мысли или критические замечания, которые у вас могут возникнуть.

Спасибо!

Ответ №1:

Я видел это раньше и надеялся, что кто-то, кто знал больше, чем я, ответил. Надеюсь, это может указать вам правильное направление

Для использования класса RectBivariateSpline вам нужны x и y в виде одномерных массивов со значениями z в виде двумерного массива (len(x), len(y))

Numpy требует, чтобы любые конкретные индексы массива были целыми числами, а не числами с плавающей точкой. Координаты x и y необходимо преобразовать в целочисленные индексы, чтобы поместить данные z в массив.

 import numpy as np
from scipy import interpolate

# Generate x and y arrays
x = np.linspace(45., 70., 70)
y = np.linspace(125/7, 59.876, 29)

# I don't know what your data looks like. Generate a list of dictionary records
data = []

for _x in x:
    for _y in y:
        data.append({'x': _x, 'y': _y, 'z': _x * _x - _y*_y/_x})

data[:100:10]
Out[4]: 
[{'x': 45.0, 'y': 17.857142857142858, 'z': 2017.9138321995465},
 {'x': 45.0, 'y': 32.86387755102041, 'z': 2000.9992344958118},
 {'x': 45.0, 'y': 47.870612244897956, 'z': 1974.075655184414},
 {'x': 45.36231884057971, 'y': 19.357816326530614, 'z': 2049.4792585649284},
 {'x': 45.36231884057971, 'y': 34.364551020408165, 'z': 2031.7068577153198},
 {'x': 45.36231884057971, 'y': 49.37128571428571, 'z': 2004.0054192006007},
 {'x': 45.72463768115942, 'y': 20.858489795918366, 'z': 2081.227345221296},
 {'x': 45.72463768115942, 'y': 35.86522448979592, 'z': 2062.610735570439},
 {'x': 45.72463768115942, 'y': 50.87195918367347, 'z': 2034.1437652565721},
 {'x': 46.08695652173913, 'y': 22.359163265306123, 'z': 2113.159976357177}]

def indexer(v):
    """ Returns a function to index into the array v by value.
        int( result   0.5 ) to avoid any not quite equal with floats.
    """
    v_lo = v.min()
    v_range = v.max()-v_lo
    n = len(v)-1

    def func( x ):
        """ Returns an index from x. """
        return int(n*(x-v_lo)/v_range 0.5)
    return func

x_index = indexer(x)  # Function to index into the x values 
y_index = indexer(y)  # Function to index into the y values

z = np.empty((len(x), len(y)))

for rec in data:
    ix_x = x_index(rec['x'])  # Map the x value to an index
    ix_y = y_index(rec['y'])  # Map the y value to an index
    z[ix_x, ix_y] = rec['z']  # Place the z value at ix_x, ix_y

inter = interpolate.RectBivariateSpline( x, y, z)

inter(45.3, 18)  # Out[6]: array([[2044.93768211]])

inter(np.arange(45., 70.), 18)
Out[7]: 
   array([[2017.8       ], [2108.95652174], [2202.10638298], [2297.25      ],
          [2394.3877551 ], [2493.52      ], [2594.64705882], [2697.76923077],
          [2802.88679245], [2910.        ], [3019.10909091], [3130.21428571],
          [3243.31578947], [3358.4137931 ], [3475.50847458], [3594.6       ],
          [3715.68852459], [3838.77419355], [3963.85714286], [4090.9375    ],
          [4220.01538462], [4351.09090909], [4484.1641791 ], [4619.23529412],
          [4756.30434783]])
  

Надеюсь, по крайней мере, некоторые указатели / идеи.

Комментарии:

1. Ваш комментарий на самом деле дал мне немного понимания, достаточного для того, чтобы заставить его работать по желанию. Я опубликую весь код через несколько часов, на случай, если кому-то это может оказаться полезным.