【问题标题】:Storing data before passing it to interpolation routine在将数据传递给插值例程之前存储数据
【发布时间】:2019-04-02 21:28:43
【问题描述】:

我将首先道歉,因为我对 python 完全陌生(这里是 Fortran 人),并且一直在学习。结果,我的知识中可能存在一些非常明显的漏洞,在阅读我目前的困境后可能会很明显。

我有一些数据需要写入文件,然后可以通过插值算法读取。目前,该插值算法可能是来自 SciPi 的 RectBivariateSpline。虽然 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 数组而不是 vanilla 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 网站上的文档并进行谷歌搜索,但到目前为止我发现的所有内容都毫无帮助。如果有人有任何关于如何使用它的好例子,那将非常有帮助。

如果您有任何建议、想法或批评,我们将不胜感激。

谢谢!

【问题讨论】:

    标签: python numpy scipy interpolation


    【解决方案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]])
    

    希望至少有一些指针/想法。

    【讨论】:

    • 您的评论确实为我提供了一些洞察力,足以让我按要求工作。我会在几个小时内发布整个代码,以防有人发现它有用。
    猜你喜欢
    • 2012-09-01
    • 2015-08-20
    • 1970-01-01
    • 2011-12-30
    • 1970-01-01
    • 2015-11-04
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多