# Line integral under surface defined by meshgrid values – Python

I need to calculate the line integral between two points (x1,y1) and (x2,y2) under a surface defined by values on a meshgrid.

I’m not exactly sure on the best tool/approach to use for this process using python.

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

```   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.
```

Any help or guidance would be appreciated.

Edit:

```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("Integral: ",I)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = plt.axes(projection="3d")

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()
```

Provided you have surface points (we can even relax the requirement of regular grid) and curve points, then basic analysis provided by `numpy` and `scipy` packages should do the trick.

First, let’s create a trial dataset for your problem.

```import numpy as np
from scipy import interpolate
```

Mainly a 3D surface:

```def f(x, y):
return x**2 + x*y + y*2 + 1

xl = np.linspace(-1.5, 1.5, 101)
X, Y = np.meshgrid(xl, xl)
Z = f(X, Y)
```

And a 2D curve:

```t = np.linspace(0, 1, 1001)
xt = t**2*np.cos(2*np.pi*t**2)
yt = t**3*np.sin(2*np.pi*t**3)
```

The complete setup looks like:

```axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.plot(xt, yt, 0)
axe.plot(xt, yt, St)
axe.view_init(elev=25, azim=-45)
```

Then we create a multidimensional linear interpolator:

```XY = np.stack([X.ravel(), Y.ravel()]).T
S = interpolate.LinearNDInterpolator(XY, Z.ravel())
```

To interpolate points from 2D curve on the 3D surface:

```xyt = np.stack([xt, yt]).T
St = S(xyt)
```

We also compute the curvilinear coordinates of the 2D curve:

```Sd = np.cumsum(np.sqrt(np.sum(np.diff(xyt, axis=0)**2, axis=1)))
```

Using curvilinear coordinates based on cumulative arc length, the integral to solve looks like:

```fig, axe = plt.subplots()
axe.plot(Sd, St[:-1])
axe.fill_between(Sd, St[:-1], alpha=0.5)
axe.grid()
```

Finally we integrate using the method of our choice, here the simplest Trapezoidal Rule from `numpy`:

```I = np.trapz(St[:-1], Sd) # 2.041770932394164
```