Issue
what I'm asking may not be possible, but I'm hoping you guys can help me out.
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.
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?
My extremely wrong attempt looks like this,
# 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))
Funny how seemingly innocuous problems can be a right pain in the a*se. :D
Edit : Here is amount of toy data,
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.]
Solution
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).
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()
Answered By - obtmind
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.