CODEWITHSUNDEEP
python-3.xnumpy
bchris
bchris

Reputation: 5

Calculate the volume of 3d plot

The data is from a measurement. The picture of the plotted data

I tried using trapz twice, but I get and error code: "ValueError: operands could not be broadcast together with shapes (1,255) (256,531)" The x has 256 points and y has 532 points, also the Z is a 2d array that has a 256 by 532 lenght. The code is below:

import numpy as np

 img=np.loadtxt('focus_x.txt')
 m=0
 m=np.max(img)
 Z=img/m 

 X=np.loadtxt("pixelx.txt",float)
 Y=np.loadtxt("pixely.txt",float)

 [X, Y] = np.meshgrid(X, Y)
 volume=np.trapz(X,np.trapz(Y,Z))

Upvotes: 0

Views: 2051

Answers (1)

gboffi
gboffi

Reputation: 25023

The docs state that trapz should be used like this

intermediate = np.trapz(Z, x)
result = np.trapz(intermediate, y)

trapz is reducing the dimensionality of its operand (by default on the last axis) using optionally a 1D array of abscissae to determine the sub intervals of integration; it is not using a mesh grid for its operation.

A complete example.

First we compute, using sympy, the integral of a simple bilinear function over a rectangular domain (0, 5) × (0, 7)

In [1]: import sympy as sp, numpy as np
In [2]: x, y = sp.symbols('x y')
In [3]: f = 1 + 2*x + y + x*y
In [4]: f.integrate((x, 0, 5)).integrate((y, 0, 7))
Out[4]: 2555/4

Now we compute the trapezoidal approximation to the integral (as it happens, the approximation is exact for a bilinear function) — we need coordinates arrays

In [5]: x, y = np.linspace(0, 5, 11), np.linspace(0, 7, 22)

(note that the sampling is different in the two directions and different from the defalt value used by trapz) — we need a mesh grid to compute the integrand and we need to compute the integrand

In [6]: X, Y = np.meshgrid(x, y)
In [7]: z = 1 + 2*X + Y + X*Y

and eventually we compute the integral

In [8]: 4*np.trapz(np.trapz(z, x), y)
Out[8]: 2555.0

Upvotes: 1

Related Questions