COMMENT, 1 day later: I now see that my solution assumes that the contours are stored as polygons. In fact, the question is about multipolygon data. For my solution to work, you will first have to extract the individual component polygons with their level, and assign them a unique integer-valued id called "id".
Two further assumptions, that are not checked, are:
- The polygons are properly nested, as they should be if they represent contours.
- The geometry field in the polygons table is called "the_geom".
Good luck! Roelof
(END OF COMMENT)
For the polygon case, I suggest creating a function using plpgsql. Here is my try:
create or replace function find_level(
point geometry, contours text default 'insert_default_here')
returns double precision as
distance_to_exterior double precision;
distance_to_interior double precision;
level_exterior double precision;
level_interior double precision;
-- If there is no polygon covering the point, the level of the outer polygon is returned:
select count(*) into point_is_covered from contours where st_contains(the_geom, point);
if point_is_covered = 0 then
return (select level from contours order by st_distance(the_geom, point) limit 1);
-- Find the id of the smallest containing polygon:
select id from contours into id_exterior where st_contains(the_geom, point) order by st_area(the_geom) limit 1;
-- If there is no smaller polygon beyond the point, the level of the containing polygon is returned:
select count(*) into interior_polygons_exist from contours where id <> id_exterior and st_within(the_geom, (select the_geom from contours where id = id_exterior));
if interior_polygons_exist = 0 then
raise debug 'No interior polygon';
return (select level from contours where id = id_exterior);
-- Find the id of the closest contained polygon:
select id from contours into id_interior where id <> id_exterior and st_within(the_geom, (select the_geom from contours
where id = id_exterior)) order by st_distance(the_geom, point) limit 1;
-- And proceed to interpolation:
select level from contours into level_exterior where id = id_exterior;
select level from contours into level_interior where id = id_interior;
select st_distance(point, st_exteriorring(the_geom)) into distance_to_exterior from contours where id = id_exterior;
select st_distance(point, the_geom) into distance_to_interior from contours where id = id_interior;
return level_exterior + distance_to_exterior / (distance_to_exterior + distance_to_interior) * (level_interior - level_exterior);
I made the assumption that your geometry field is called the_geom. You may have to substitute this value with the actual name in your table.
To apply the function:
select *, find_level(the_geom, 'name_of_contour_table')
For the contour lines, it's harder. If possible, do an intermediate step to convert the contour lines to polygons.