title: Line integral under surface defined by meshgrid values - Python
body: <p>I need to calculate the line integral between two points (x1,y1) and (x2,y2) under a surface defined by values on a meshgrid.</p>
  <p><a href="https://i.stack.imgur.com/nX7PO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/nX7PO.png" alt="enter image description here" /></a></p>
  <p>I'm not exactly sure on the best tool/approach to use for this process using python.</p>
  <p>As I do not have a function which represents the surface, instead values at points on a evenly spaaced meshgrid I am assuming I will need to use one of the following methods</p>
  <pre><code>   trapz         -- Use trapezoidal rule to compute integral from samples.
     cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
     simps         -- Use Simpson's rule to compute integral from samples.
     romb          -- Use Romberg Integration to compute integral from
                      (2**k + 1) evenly-spaced samples.
  </code></pre>
  <p>Any help or guidance would be appreciated.</p>
  <p>Edit:</p>
  <pre><code>import numpy as np
  from scipy import interpolate
  
  def f(x, y):
      return x**2 + x*y + y*2 + 1
  
  xl = np.linspace(-1.5, 1.5, 101,endpoint = True)
  X, Y = np.meshgrid(xl, xl)
  Z = f(X, Y)
  
  #And a 2D Line:
  arr_2D = np.linspace(start=[-1, 1.2], stop=[0, 1.5], num=101,endpoint = 
  True) #Creates a 2D line between these two points
  
  #Then we create a multidimensional linear interpolator:
  XY = np.stack([X.ravel(), Y.ravel()]).T
  S = interpolate.LinearNDInterpolator(XY, Z.ravel())
  print(S)
  
  #To interpolate points from 2D curve on the 3D surface:
  St = S(arr_2D)
  
  #We also compute the curvilinear coordinates of the 2D curve:
  
  #Using curvilinear coordinates based on cumulative arc length, the integral to solve looks like:
  Sd = np.cumsum(np.sqrt(np.sum(np.diff(arr_2D, axis=0)**2, axis=1)))
  print(Sd)
  
  I = np.trapz(St[:-1], Sd) # 2.041770932394164
  print(&quot;Integral: &quot;,I)
  
  import matplotlib.pyplot as plt
  from mpl_toolkits.mplot3d import Axes3D
  
  fig = plt.figure()
  ax = plt.axes(projection=&quot;3d&quot;)
  
  x_line = np.linspace(start=[-1], stop=[1.5], num=100,endpoint = True)
  y_line = np.linspace(start=[-1.2], stop=[1.5], num=100,endpoint = True)
  
  ax.plot3D(x_line, y_line, 'red')  #Line which represents integral
  ax.plot_wireframe(X, Y, Z, color='green') #Represents the surface
  ax.set_xlabel('x')
  ax.set_ylabel('y')
  ax.set_zlabel('Time')
  
  plt.show()
  
  fig = plt.figure()
  ax = plt.axes()
  ax.fill_between(Sd, St)
  ax.set_xlabel('x')
  ax.set_ylabel('Z')
  plt.show()
  </code></pre>