1

I would like to render contours using Mapnik. The problem is that some contours are very long which makes rendering too slow and so I decided to spling long contours to multiple smaller segments. Input table (generated with gdal_contour):

gis=> \d contour;
                                            Table "public.contour"
    Column    |            Type             | Collation | Nullable |                 Default                  
--------------+-----------------------------+-----------+----------+------------------------------------------
 ogc_fid      | integer                     |           | not null | nextval('contour_ogc_fid_seq'::regclass)
 id           | numeric(8,0)                |           |          | 
 height       | numeric(12,3)               |           |          | 
 wkb_geometry | geometry(LineString,900914) |           |          | 
Indexes:
    "contour_pkey" PRIMARY KEY, btree (ogc_fid)
    "contour_wkb_geometry_geom_idx" gist (wkb_geometry)

I am no GIS or PostGIS proffesional but I tried my best and came up with following query:

create table contour_split as select
  ogc_fid,
  height,
  (select st_makeline((x).geom) from (select st_dumppoints(wkb_geometry) as x offset round(start) limit round(start + seg_len) - round(start)) as sub3) as geom
from
  (
    select
      ogc_fid,
      wkb_geometry,
      height,
      n_pts / ceil(n_pts / 1000.0) as seg_len,
      generate_series(
        0,
        n_pts - n_pts / ceil(n_pts / 1000.0),
        n_pts / ceil(n_pts / 1000.0)
      ) as start
    from
      (select *, st_numpoints(wkb_geometry) as n_pts from contour) as sub1
  ) as sub2;

The problem is that this query is very slow. There are over 1.6M rows in the contours table and the query is running for several hours already.

I've also tried to do it with intermediate tables but the speed (of the last query) is still slow:

alter table contour add column np int;

update contour set np = st_numpoints(wkb_geometry);

create table contour_tmp as select
    ogc_fid,
      generate_series(
        0,
        np - np / ceil(np / 1000.0),
        np / ceil(np / 1000.0)
      ) as start
    from contour;

create table contour_split as
select
  ogc_fid
  height,
  (
    select
      st_makeline((x).geom)
    from
      (
        select
          st_dumppoints(wkb_geometry) as x offset round(start)
        limit
          round(start + np / ceil(np / 1000.0)) - round(start)
      ) as ddd
  )
from
  contour
  join contour_tmp using (ogc_fid);

Can somebody see there something "obvious" what could improve the speed of the query? Is there a better approach to achieve even splitting linestrings to not to be longer than 1000 points? ST_Subdivide returns uneven results - some segments are very short. Or to make Mapnik fast even with long contour linestrings?

Version information: POSTGIS="2.5.1 r17027" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" GDAL="GDAL 2.4.0, released 2018/12/14" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" RASTER

1

I've created a simple program that does what I need in less than a hour:

const { Client } = require('pg');
const wkx = require('wkx');
const QueryStream = require('pg-query-stream');

const client = new Client();
const client2 = new Client();

let num = 0;

const fn = async () => {
  await Promise.all([
    client.connect(),
    client2.connect(),
  ]);

  const qs = new QueryStream('SELECT ogc_fid, height, st_asbinary(wkb_geometry) AS wkb_geometry FROM contour');
  const stream = client.query(qs);

  for await (const row of stream) {
    num++;
    if (num % 1000 === 0) {
      console.log('ROW', num);
    }
    const gj = wkx.Geometry.parse(row.wkb_geometry).toGeoJSON();
    // console.log(gj);
    const len = gj.coordinates.length;
    const n = Math.ceil(len / 1000);
    const size = len / n;
    let from = 0;

    for (let i = 0; i < n; i++) {
      const rSize = Math.round(size);
      const sliceCoords = gj.coordinates.slice(from, from + rSize);
      const sliceGeom = new wkx.LineString(sliceCoords.map(([x, y]) => (new wkx.Point(x, y)))).toWkb();
      await client2.query(
        'INSERT INTO contour_split (cid, height, geom) VALUES ($1, $2, ST_GeomFromWKB($3, 900914))',
        [
          row.ogc_fid,
          row.height,
          sliceGeom,
        ],
      );

      from += rSize - 1;
    }
  }

  client.end();
  client2.end();
};

fn().catch(err => {
  console.log(err);
});
0

The PostGIS ST_Subdivide function is the command you are looking for. It will allow you to set a vertex limit, after which the feature is "subdivided" / cut in two, so no feature will have more than X vertices.

  • Please see my note in the question: ST_Subdivide returns uneven results - some segments are very short – Martin Ždila Jun 24 at 14:53
  • That is not really relevant, the number of vertices should be the same (at least, I assume there is no serious bug in the function...). If some segments are very short, they are likely related to contour segments with many vertices, e.g. contours in a DEM that hasn't been cleaned for forests, might yield very detailed contours with many bends and curves for individual canopies of trees. – Marco_B Jun 24 at 15:04
  • Another option is that gdal_contour has an issue and returns contours with uneven number of vertices (seems unlikely as well though). – Marco_B Jun 24 at 15:05
  • I mean that ST_Subdivide divides linestring to segments of with "random" number of vertices. Example: pastebin.com/tLTb6Dmm – Martin Ždila Jun 24 at 15:21
  • Ah, to be honest, I never really checked this... I see two possibilities for this to happen: 1) PostGIS uses a in the Help poorly documented routine that cuts each line exactly in half - BASED ON LENGTH! - until each segment falls below the thresshold in an attempt to satisfy a secondary condition to not create segments differing to much in length (as you can imagine, otherwise some segments might be much longer than others in terms of actual length based on the detail of the contour in a certain regions of the extent of the dataset). – Marco_B Jun 24 at 15:32

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.