6

Using PostGIS, what is the best method to interpolate the z value for points (the orange dot in the image below) located between each contour layer (blue lines and red lines in the image below)? So for the blue contours it's going to be ~ 917 and for the red contours it'll be ~ 819.

example-map

The blue contours are a 2D Multipolygon and the red contours are a 2D Multilinestring (so for both, the contour value is an attribute field. If necessary I could update the geometries to have a z value).

Thanks!

1

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:

  1. The polygons are properly nested, as they should be if they represent contours.
  2. 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
$body$
    declare
        point_is_covered integer;
        interior_polygons_exist integer;
        id_exterior integer;
        id_interior integer;
        distance_to_exterior double precision;
        distance_to_interior double precision;
        level_exterior double precision;
        level_interior double precision;
    begin
        -- 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);
        end if;

        -- 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);
        end if;

        -- 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);
    end;
$body$
language plpgsql;

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.

0

There might be a way to create rasters from these multipolygons if they were to include z-values and then access individual points on them using ST_PixelAsPoint

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.