Python function to find the numeric volume integral?

Goal

I would like to compute the 3d volume integral of a numeric scalar field.

Code

For this post, I will use an example of which the integral can be exactly computed. I have therefore chosen the following function:

In python, I define the function, and a set of points in 3D, and then generate the discrete values at theses points:

import numpy as np


# Make data.
def function(x,y,z):
    return x**y**z

N = 5
grid = np.meshgrid(
    np.linspace(0, 1, N),
    np.linspace(0, 1, N),
    np.linspace(0, 1, N)
    )

points = np.vstack(list(map(np.ravel, grid))).T

x = points[:,0]
y = points[:,1]
z = points[:,2]

values = [function(points[i,0],points[i,1], points[i,2]) for i in range(len(points))]

Question

How can I find the integral, if I don't know the underlying function, i.e. if I only have the coordinates (x,y,z) and the values?

from Python function to find the numeric volume integral?