Computing volume of concave polyhedron

Question

I have a circular grid with points uniformly distributed throughout it. See this:

Each point has some nonnegative height assigned to it (i.e. height can be 0 on up).

I'm trying to accurately compute the volume of this...polyhedron(?)

My work so far has been looking clusters of 4 points (i.e. little squares inside the grid made up of 4 points), and trying to calculate the volume of the polyhedron made up of those 4 points and the 4 points with the height values. Think this:

The trouble with this method is that the polyhedrons aren't necessarily convex, and that means I can't easily use an algorithm to break up each polyhedron into tetrahedron and sum the their volumes.

Any thoughts?

Would it be better to apply some sort of algorithm to the entire circular grid which could break it into disjoint, convex polyhedrons; then use another algorithm to compute the volume of each one via breaking it into tetrahedrons?

I apologize for what may be a noobish question, and for not having much in the way of work to show what I have attempted. Geometry has never been a strength of mine, and I'm not really sure the best way to start approaching this problem.

edit/update Thanks for all the quick answers / suggestions! This community always amazes me. Sorry for my delay in responding. Went out of town today to celebrate my niece's birthday. I'll do my best to address some of the questions left in the comments and help to clarify the specifics of my problem. Apologies in advance for my ignorance about the geometry at play here. I'm a statistician by schooling, and never did great with geometry.

1. The circle is a region ~38.4 square km area centered around a GHCN weather station. The points inside the circle are evenly spread, and they represent locations where precipitation values were interpolated using ordinary Kriging, the height associated with each point (shown in the second picture) is the precipitation value in 10mm increments. The goal is measuring the volume of precipitation over that circular region.

My initial, and what I assumed was going to be error prone, method involved dividing the total area of the circular region by the number of points in the grid, taking that value as the area for the region centered around each point. The height at this point was multiplied the the area, and all of the resulting volumes were summed. (Henning, I think this is similar to what you proposed in the comments.) However, I figured there was a better method out there that I might be able to implement in R or matlab.

2. The object in the second picture...I got the idea to look at the squares of 4 points in the grid, the heights at each point, and using these coordinates to define a the vertices of a polyhedron. If I could compute the volume of each of these, then I can write a function to go over that grid and compute the entire volume of precipitation more accurately (I think) than the first method. The second picture was just a quick 3D plot I made in R to help visualize what one such polyhedron could look like. That was when I posted here.

3. Studiosus- your solution seems rather promising because I think I could easily define the edges, faces, and vertices. As for reorienting them and doing the integration...do you have any resources you could point me to where I could read up on how that all works? I'm not sure I could dive in an code a function to do all of that quite yet.

4. Heimdall and JeanMarie- this solution seems perfect, and quite easy to implement computationally. To make sure I didn't gloss over anything, I would sum the 4 heights for each 2x2 square, divide by 4, and multiply by the area of the square (as the formula given uses a unit square). Do this for each square and sum. Right? If so, this seems like my best bet.
   Follow up question: some boundary cases crop up where I have 3 points making a triangular base. The volume for these would be computed by slicing the object as a triangular prism, with height equal to the shortest of the 3 heights, and the remaining object would be a tetrahedron. I can compute both those volumes easily and add them together. Does this seem like a good way to handle those cases?

5. Finally, am I gaining a lot of accuracy by employing the method given by JeanMarie and Heimdall over my initial method? Would there be a reason to pursue the solution proposed by Studiosus instead of JeanMarie and Heimdall given the context of my problem?

Thank y'all so much for the advice and help so far!

Answer 1

I approve @Heimdall.

Let me explain why such a simple formula is valuable.

You may know that one of the best local approximation of a surface built on a square grid like yours is through bilinear interpolation : let us consider the simple case of points with coord. $(0,0,h_{00}),(1,0,h_{10}),(1,1,h_{11}),(0,1,h_{01})$. Then the quadric surface passing through these 4 points (in general a PH, abbreviation for Parabolic Hyperboloid, see figure below) has the following equation:

$$z=f(x,y)=(1-x)(1-y)h_{00}+x(1-y)h_{10}+(1-x)yh_{01}+xyh_{11}$$

(check for example that, when $x=1$ and $y=0$, the only nonzero "surviving" term is the height $h_{10}$ precisely associated with point $(1,0,h_{10})$).

Now, the volume under this surface is easily computed to be $$\int_{x=0}^1\int_{y=0}^1 \ f(x,y)dxdy=\dfrac{h_{00}+h_{10}+h_{01}+h_{11}}{4}$$

Now, adapting the above formula to the grid's steps, is easy and amounts, as it has been said, to the multiplication of the above formula by the area of the grid elementary square. Dimensionaly speaking : meters by square meters make cubic meters, the right unit for a volume; everything is for the best...