[postgis-users] Z Interpolation of Point in Triangle

[Steffen Macke]

Hello All,

the PL/SQL function below allows a "poor man's" TIN implemenation,
with a triangle represented by a polygon (3 dimensions). The function
returns the interpolated z value
if it is fed with such a triangle polygon and a point that lies inside
the polygon.

It would be nice to have it included in the PostGIS distribution.

An example how to use the function is included at the end.

Regards,

Steffen

create or replace function interpolate_from_tin(location geometry, tin
geometry) returns float as
$body$
declare
  az float8;
  bz float8;
  cz float8;
  a1 float8;
  a2 float8;
  b1 float8;
  b2 float8;
  c1 float8;
  c2 float8;
  result float8;
  ring geometry;
  tinpoint geometry;
  a geometry;
  b geometry;
  c geometry;
begin
  if numgeometries(tin) <> 1 then
    raise exception 'Expecting 1 geometry in TIN polygon, got %.',
numgeometries(tin);
  end if;
  ring := exteriorring(geometryn(tin,1));
  if numpoints(ring) <> 4 then
    raise exception 'Expecting 3 points in TIN polygon, got %.',
numpoints(ring);
  end if;
  if ndims(ring) < 3 then
    raise exception 'Expecting 3 dimensions in TIN polygon, got %.',
ndims(ring);
  end if;
  a := pointn(ring, 1);
  b := pointn(ring, 2);
  c := pointn(ring, 3);
  a1 := x(location)-x(a);
  a2 := y(location)-y(a);
  b1 := x(location)-x(b);
  b2 := y(location)-y(b);
  c1 := x(location)-x(c);
  c2 := y(location)-y(c);
  az := z(a);
  bz := z(b);
  cz := z(c);

  result := ((a1*b2*cz)-(a1*bz*c2)+(a2*bz*c1)-(a2*b1*cz)+(az*b1*c2)-(az*b2*c1))/((a1*b2)-(a1*c2)+(a2*c1)-(a2*b1)+(b1*c2)-(b2*c1));

  return result;
end;
$body$
language plpgsql immutable;

comment on function interpolate_from_tin(geometry, geometry) is
 'Linear interpolation of a given points z value from a triangle polygon
 @param location is the point geometry.
 @param tin is the triangle polygon (3D).
 @returns the interpolated z value at point location.';

select interpolate_from_tin('POINT(1 1)', 'MULTIPOLYGON(((0 0 0,1 2
1,2 0 0,0 0 0)))'),
 interpolate_from_tin('POINT(1 2)', 'MULTIPOLYGON(((0 0 0,1 2 1,2 0
0,0 0 0)))'),
 interpolate_from_tin('POINT(2 0)', 'MULTIPOLYGON(((0 0 0,1 2 1,2 0
0,0 0 0)))');