Split long lines in PostGIS

Question

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

Answer 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);
});

Answer 2

Another approach may be to use ST_Segmentize(geom, length), which allows to split Linestrings in segments not longer than length.

Answer 3

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.