Surface plot with multiple polynomial fits
<p>what I'm asking may not be possible, but I'm hoping you guys can help me out.</p>

<p>So I have two 2D arrays, f1(x1) = y1, and f2(x2) = y2. I want to make a surface plot of the ratio of these, so the z dimension is (f2(x2)/f1(x1)). Unfortunately I'm coming up against a wall from whatever direction I approach the problem.</p>

<p>My main problem is that the ranges of each array is different, x1 goes from 300 to 50,000, and x2 goes from 300, 1200. Now I'm happy with assuming that f2(x2) = f2(1200) for all x2 > 1200. But this bound means it's impossible for me to fit a polynomial this data in any realistic way (my first set of data is nicely reproduced by a 5th order polynomial, and my second set of data is best fit with a 1st order polynomial). Is there an alternative way that I can fit a function to (x2,y2) so it takes the boundary values for all points outside of the boundaries?</p>

<p>My extremely wrong attempt looks like this,</p>

<pre><code># x1 is an array from 300 to 50,000 in steps of 50
# x2 is an array from 300 to 1,150 in steps of 50

f1_fit = np.poly1d(np.polyfit(x1, y1, 5))
f2_fit = np.poly1d(np.polyfit(x2, y2, 1))

X, Y = np.meshgrid(x1, x2)
Z = (f2_fit(x2) / f1_fit(x1))
</code></pre>

<p>Funny how seemingly innocuous problems can be a right pain in the a*se. :D</p>

<p>Edit : Here is amount of toy data,</p>

<pre><code>x1 = x2 = [  300.   350.   400.   449.   499.   548.   598.   648.   698.   748.
             798.   848.   897.   947.   997.  1047.  1097.  1147.  1196.  1246.
             1296.  1346.  1396.  1446.  1496.  1546.  1595.  1645.  1695.  1745.]

y1 = [  351.   413.   476.   561.   620.   678.   734.   789.   841.   891.
        938.   982.  1023.  1062.  1099.  1133.  1165.  1195.  1223.  1250.
        1274.  1298.  1320.  1340.  1360.  1378.  1395.  1411.  1426.  1441.]

y2 = [ 80.  75.  70.  65.  62.  58.  58.  52.  48.  46.  44.  41.  38.  35.  32.
       32.  29.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.]
</code></pre>

<p>So I managed to solve my problem. I did some preproccessing on my data as explained above, by setting x1 = x2, and thus extrapolated the edge values for f(x2).</p>

<pre><code>import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

data1 = np.loadtxt('data1.dat')
data2 = np.loadtxt('data2.dat')

X = []
Y = []
Z = []

for i in data1:
    for j in data2:
        X.append(i[0])
        Y.append(j[0])
        Z.append(i[1]/j[1])

x_mesh,y_mesh, = np.meshgrid(np.linspace(300,50000,200), np.linspace(300,50000,200))

z_mesh = interp.griddata((X,Y),Z,(x_mesh,y_mesh),method='linear')

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(x_mesh,y_mesh,z_mesh,cstride=10,rstride=10,cmap='OrRd_r')
plt.show()
</code></pre>

python matplotlib curve-fitting surface