Create contour line from point table in PostGIS

Question

I have a point table containing values ranging from 0 to 100 and I would like to draw a contour line for the value 70. (See image below)

I know how to create contour lines in QGIS for raster data but now I would like to do it with vector data (my point table) using PostGIS so that it can be integrated in a python script.

Any idea on how to start?

EDIT: This question is different from Generating contour lines and filled contours from points as it doesn't involve raster data (neither in the input nor in the output)

Answer 1

I think what you want is something like gdal_contour.

gdal_contour: Builds vector contour lines from a raster elevation model.


Usage: gdal_contour [-b <band>] [-a <attribute_name>] [-3d] [-inodata]
                    [-snodata n] [-i <interval>]
                    [-f <formatname>] [[-dsco NAME=VALUE] ...] [[-lco NAME=VALUE] ...]
                    [-off <offset>] [-fl <level> <level>...]
                    [-nln <outlayername>]
                    <src_filename> <dst_filename>

But PostGIS doesn't seem to implement anything like that from GDAL.

GDAL expects a input raster elevation model (DEM). PostGIS has a few DEM functions but I don't believe they'll help you here.

Answer 2

To create contour in PostGIS I usually use ST_ConvexHull or ST_ConcaveHull

Usage of it is simple as that:

SELECT
  st_convexhull(st_collect(geom)),
               pointsId
             FROM point_table 
             GROUP BY pointId;

But to create polygon like at your picture is not really easy.

You can try some other functions using this cheatsheet.

Answer 3

I am trying to do something very similar, I found an item at [fast contour lines using plain postgis1. it creates a contour map from your list of points. I get the impression it expects rectangle input.

My problem with this solution was the calling method, it uses array of data, a concept which is outside of most sql grammers.

My data weather point data, where geom the point location, speed a scalar value denoting the current wind speed and reading_time the value data was collected.

My calling code is

drop table foo;
create table foo(
    id serial,
    geom geometry ,
    break numeric 

);
insert into foo (geom,break )  
select 
    (contour_lines(  array_agg(geom),array_agg(speed::numeric), 
    ARRAY[ 0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5,14.0,14.5,15.0,15.5,16.0,16.5] )).geom ,  
    (contour_lines(  array_agg(geom),array_agg(speed::numeric), 
    ARRAY[ 0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5,14.0,14.5,15.0,15.5,16.0,16.5] )).break
from windproc where reading_time='2018-08-28 07:00:00+01'  and geom is not null and speed is not null;

Where array_agg creates an array of the object supply as an argument. The array construct creates a list of step/break values which is used a banding values for the generated contour lines.

The code for contour_lines is provided by the link. it expects a numeric value for scalar value (in my case speed) hence the cast.

The reason there is what appears to be a double called to contour_lines is that you need to get the break and geom answers provide by contour_lines function, but since array_agg is an aggrate function it can only appear in some parts of sql construct.

The result from this will give you a table called foo with the desired information in the geom column