Thread Subject:

how to calculate area of interpolated surface

I want to compute the area of a interpolated surface. I used the "interp2" method so I got no function relation between the (x,y) and z points.

Does any of you have a suggestion how to do this?

Thanks.

"Tomas " <tomasnic@gmail.com> wrote in message <gom2b9$m0n$1@fred.mathworks.com>...
> I want to compute the area of a interpolated surface. I used the "interp2" method so I got no function relation between the (x,y) and z points.
>
> Does any of you have a suggestion how to do this?
>
> Thanks.

Since you used interp2, I'll assume that you
have a lattice of points on a regular grid, and
therefore a function evaluated at the lattice of
points f(x,y).

Now, to compute the area of that surface, there
are several simple solutions. The first is to use
a rectangle rule. So assume we have a function
F as an nxm array, evaluated on the nxm lattice
of points. Assume the spacing in x and y is known
and fixed, at dx and dy respectively.

Assume the surface is defined as a series of
rectangles, at the height of one of the corners
of each square in the lattice. Then the surface
area is given as

  sarea1 = (n-1)*(m-1)*dx*dy

This is the rectangle rule integration. As you can
see, it misses one row and column entirely. It
very much depends on how you will define the
surface though. If instead, you wish to define
the surface as a series of rectangles, each
centered exactly on each lattice point, then the
surface area is simply

  sarea2 = n*m*dx*dy;

dx = 3;
dy = 2;
n = 10;
m = 20;
F = rand(n,m);
surf(F)

Most likely, both of these approximations are
not adequate for your purposes. So a better
solution might be to form a set of triangles.
Triangulate the lattice, breaking each square
into a pair of triangles. Find the area of each
triangle, then sum over the areas. It is doable,
and perhaps before I get done writing I'll do
exactly that.

sarea1 = (n-1)*(m-1)*dx*dy
sarea1 =
        1026

sarea2 = n*m*dx*dy
sarea2 =
        1200

(Ok, I did it just now, but the method I actually
used was not easy to explain, and would just
have confused things.) Here is the area of the
triangles. As you should expect, it is in the
middle of the others.

satri
satri =
       1170.7

Here is another estimate, formed by approximating
the surface area by integrating the length of the
piecewise linear line segment along one of the
two dimensions.

trapz(sum(dx*sqrt(1 + diff(F,[],1).^2),1),2)*dy
ans =
         1099

I'll choose to do the same in the other dimension
too, just to see if it made a difference.

trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx
ans =
       1101.9

Happily the difference is slight. The point is,
your solution must depend upon how you define
the surface. There are other, higher order
approximations one might choose for all of this.

I'll just stop here.

John

Hi John,

thanks for your answer, but the problem is that I got no function z = F(x,y) at the lattice points.

Maybe it will help if I describe what I want to do. I want to compute the area of a small region of a topographical map. Maybe I can fit some kind of polynomial function to points on the contour lines using interpolation or a overestimation method. When I got a function I will try the method you described.

Tomas

"Tomas " <tomasnic@gmail.com> wrote in message <goosuj$nr0$1@fred.mathworks.com>...
> Hi John,
>
> thanks for your answer, but the problem is that I got no function z = F(x,y) at the lattice points.
>
> Maybe it will help if I describe what I want to do. I want to compute the area of a small region of a topographical map. Maybe I can fit some kind of polynomial function to points on the contour lines using interpolation or a overestimation method. When I got a function I will try the method you described.
>
> Tomas

Then what do you have? All that I assumed is that
you have a list of values, defined at each point on
a regular lattice. If you used interp2, then you MUST
have that much.

Do you have points only on the contour lines? If
so, then how do you claim to have used interp2?
This would be impossible.

So tell me what exactly you do have? Be precise
when you ask a question, else all you do is waste
the time of someone who might try to answer the
wrong question (as apparently did I.)

John

"Tomas " <tomasnic@gmail.com> wrote in message <goosuj$nr0$1@fred.mathworks.com>...
> .......
> I want to compute the area of a small region of a topographical map.
> .......

  Tomas, I am guessing you want to find the approximate surface area of a region on your map, as opposed to the area of that region as projected onto a horizontal plane. That is, you want the "slant" area which takes into account the gradients present in the region. And you wish to use the output of 'interp2' for this purpose.

  It should be pointed out that in the case of mountainous terrains such a surface figure would be highly dependent on the fineness of resolution of your data. Down at a fine resolution the surface area may be vastly greater than at a course resolution because of various surface irregularities. (As a one-time backpacker I can testify to the truth of that statement!)

  In calculus for an infinitesimal rectangular region dx by dy, its three-dimensional surface area is given by

 sqrt(1 + (dz/dx)^2+(dz/dy)^2)*dx*dy

You can therefore use finite differences for approximating this expression for each small rectangular region in the mesh of points XI,YI,ZI to get a first order approximation to the desired surface area. Assuming your data has a reasonably high resolution, that should be adequate for your purposes.

Roger Stafford

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gom6nt$8m1$1@fred.mathworks.com>...
> .......
> trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx
> ans = 1101.9
> .......
> John

  Note to John: I see you have also done a surface calculation above but I think that gets an answer that is too small. Shouldn't it be changed to something like this:

 trapz(sum(dy*sqrt(1 + diff(F,[],1).^2 + diff(F,[],2).^2),2),1)*dx

to get an accurate multiplication by the true secant of the angle between the vertical and the surface normal? Or is my head not working properly this morning?

Roger Stafford

"John D'Errico" <woodchips@rochester.rr.com> wrote in message
> Then what do you have? All that I assumed is that
> you have a list of values, defined at each point on
> a regular lattice. If you used interp2, then you MUST
> have that much.
>
> Do you have points only on the contour lines? If
> so, then how do you claim to have used interp2?
> This would be impossible.
>
> So tell me what exactly you do have? Be precise
> when you ask a question, else all you do is waste
> the time of someone who might try to answer the
> wrong question (as apparently did I.)
>
> John

I?m very sorry wasting your time, John. I realize now that I can?t use interp2 because i don?t have function values for an equidistant grid. I promise to think more about my questions before posting later.

Roger, thanks for your answer. I now have an idea what to do.

Tomas

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gop9sd$3ib$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gom6nt$8m1$1@fred.mathworks.com>...
> > .......
> > trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx
> > ans = 1101.9
> > .......
> > John
>
> Note to John: I see you have also done a surface calculation above but I think that gets an answer that is too small. Shouldn't it be changed to something like this:
>
> trapz(sum(dy*sqrt(1 + diff(F,[],1).^2 + diff(F,[],2).^2),2),1)*dx

Not quite, though I was sloppy. Best is probably this:

[X,Y] = meshgrid(0:.1:1,0:.2:3);
dx = 0.1;
dy = 0.2;
[Fx,Fy] = gradient(F,dx,dy);
trapz(trapz(sqrt(1+Fx.^2 + Fy.^2)))*dx*dy
ans =
       46.689

Which correctly integrates the area, using gradient
to compute the partials.

This generates a surface area that is quite consistent
with an area computed from the sum of the triangulated
area as 46.811.

I don't know how to solve the problem as posed by the
OP, since we still don't know what he has.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gopjc6$4i0$1@fred.mathworks.com>...
> I don't know how to solve the problem as posed by the
> OP, since we still don't know what he has.

I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.

Tomas

"Tomas " <tomasnic@gmail.com> wrote in message <gotbpa$prr$1@fred.mathworks.com>...
> I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.
> .......

  I'm just speculating here, but suppose you have a way of obtaining the varying "orthogonal" distance between adjacent contour lines. If h is the fixed height difference represented between the contours and d the orthogonal horizontal distance represented between them, then the slant distance across between the contours is just sqrt(d^2+h^2). Integrating this varying slant distance along the space between the contours gives the surface area it represents, as opposed to the horizontally projected area that appears on your map between the contours. Repeating this for all spacings between contours for a region would give you its total surface area.

  Of course this is only approximate, but then being restricted to discrete contour lines is necessarily an inexact procedure. However I am not sure just how you would go about determining this variable quantity d above starting with just a contour map. If the contours are close together they would ordinarily be nearly parallel but widely-spaced contours can sometimes be decidedly non-parallel. Fortunately this latter is just the case where the needed correction over projected area is the smallest.

  Well this is just a wild idea for you to play with.

Roger Stafford

"Tomas " <tomasnic@gmail.com> wrote in message <gotbpa$prr$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gopjc6$4i0$1@fred.mathworks.com>...
> > I don't know how to solve the problem as posed by the
> > OP, since we still don't know what he has.
>
> I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.
>
> Tomas

As Roger pointed out, we need to be able to evalute the gradient of the topographic map at any point in order to compute the area. Therefore I would suggest to approximate your surface map by the spline surface - it is twice derivable - so it is nice to integrate numerically its gradient. That might be violate the property of the map, because land might have vertical fall (infinity gradient), up to you...

To find the spline surface, you could create a regular grid. Next find find out what to fill on the grid in order to have the contour line data. This is a linear problem : you can solve it by building the matrix (using interp2) and using "\" operator.

Now there is a little bit of work to extract the gradient at a given point. You can either do your own calculation from the spline data values. You might want to know Matlab use a not-a-knot condition for interpolation to calculate for points at the outer ring. You could also use 1D spline functions applied to 1D grid crossing the point of interest (that gives you a polynomial coefficient on all braket that you can take the derivative). I believe the spline toolbox have some nice feature to do this kind of calculation.

The last step: The calculation of the area by using function quad2d. You have to program a function that returns 0 outside the region, and in infinitisimal surface area (see Roger's above formula) inside the region - now you have everything to calculate it. You could also do a brute force sum of infinitisimal surface areas with point that distributed regularly inside the region. You might need a tool to tell when a point is inside the polygonal contour line. I believe there are several such functions on FEX. Be aware about the regions with hole inside it that need some care in your implementation.

Good luck

Bruno