n50_path.sql

-- Shortest path routing function for N50 vegsti data

CREATE OR REPLACE FUNCTION path(
  IN x1 double precision,
  IN y1 double precision,
  IN x2 double precision,
  IN y2 double precision,
  IN path_buffer double precision DEFAULT 2000,
  IN point_buffer double precision DEFAULT 10,
  IN targets integer DEFAULT 1,
  IN srid_in integer DEFAULT 4326,
  IN srid_db integer DEFAULT 25833,
  IN bbox double precision ARRAY[4] DEFAULT '{}',
  OUT path_id integer,
  OUT cost double precision,
  OUT geom geometry
) RETURNS SETOF record AS

$BODY$
DECLARE
  sql       text;
  point1    geometry;
  point2    geometry;
  bounds    text;
  rec       record;
  source    record;
  target    record;
  prec1     double precision;
  prec2     double precision;

BEGIN
  -- Find the closest edge (source) near the start (x1, y1)
  point1 := ST_Transform(ST_SetSRID(ST_MakePoint(x1, y1), srid_in), srid_db);

  EXECUTE 'SELECT
      ogc_fid AS id,
      ST_LineLocatePoint(geometri, $1) AS prec
    FROM n50.n50_vegsti
    WHERE geometri && ST_Buffer($1, ' || point_buffer || ')
    ORDER BY ST_Distance(geometri, $1)
    LIMIT 1'
  INTO source
  USING point1;

  -- Find the closest edge (target) near the end (x2, y2)
  point2 := ST_Transform(ST_SetSRID(ST_MakePoint(x2, y2), srid_in), srid_db);

  EXECUTE 'SELECT
      ogc_fid AS id,
      ST_LineLocatePoint(geometri, $1) AS prec
    FROM n50.n50_vegsti
    WHERE geometri && ST_Buffer($1, ' || point_buffer || ')
    ORDER BY ST_Distance(geometri, $1)
    LIMIT 1'
  INTO target
  USING point2;

  RAISE NOTICE '[ROUTER] source.id=% source.prec=%', source.id, source.prec;
  RAISE NOTICE '[ROUTER] target.id=% target.prec=%', target.id, target.prec;

  -- Return if we could not find source or target edges
  IF source.id IS null OR target.id IS null THEN
    RAISE NOTICE '[ROUTER] source or target is NULL; returning';
    RETURN;
  END IF;

  -- Since there are millions of edges in the database we need to limit the
  -- shortes path search to only the subset of edges that can reasonably match.
  IF array_length(bbox, 1) != 0 THEN
    RAISE NOTICE '[ROUTER] creating bounds based on bbox';
    bounds := 'ST_Transform(ST_MakeEnvelope(
        ' || bbox[1] || ', ' || bbox[2] || ',
        ' || bbox[3] || ', ' || bbox[4] || ',
        ' || srid_in || '
    ), ' || srid_db || ')';
  ELSE
    RAISE NOTICE '[ROUTER] creating bounds based on path_buffer';
    bounds := 'ST_Buffer(ST_Transform(ST_MakeEnvelope(
      ' || x1 || ', ' || y1 || ', ' || x2 || ', ' || y2 || ', ' || srid_in || '
    ), ' || srid_db || '), ' || path_buffer || ')';
  END IF;

  -- Construct a query to find all edges that should be routed for shortest
  -- path. This is not the shortest path routing but a set of edges that the
  -- router will work on.
  sql := 'SELECT
    ogc_fid as id,
    source::int,
    target::int,
    cost::float
  FROM n50.n50_vegsti
  WHERE geometri && ' || bounds || '

  UNION SELECT
    ogc_fid AS id,
    source::int,
    -888::int AS target,
    cost::float * ' || source.prec || ' AS cost
  FROM n50.n50_vegsti
  WHERE ogc_fid = ' || source.id || '

  UNION SELECT
    ogc_fid AS id,
    -888::int AS source,
    target::int,
    cost::float * (1 - ' || source.prec || ') AS cost
  FROM n50.n50_vegsti
  WHERE ogc_fid = ' || source.id || '

  UNION SELECT
    ogc_fid AS id,
    source::int,
    -999::int AS target,
    cost::float * ' || target.prec || ' AS cost
  FROM n50.n50_vegsti
  WHERE ogc_fid = ' || target.id || '

  UNION SELECT
    ogc_fid AS id,
    -999::int AS source,
    target::int,
    cost::float * (1 - ' || target.prec || ') AS cost
  FROM n50.n50_vegsti
  WHERE ogc_fid = ' || target.id;

  -- If the edge is long or the route is short we could end up with a shortest
  -- route consisting of a single edge or a subset of it. We handle that by
  -- inserting dummy edges into the pool of edges.

  IF (source.id = target.id) THEN
    RAISE NOTICE '[ROUTER] source.id = target.id';

    sql := sql || 'UNION SELECT
      ogc_fid AS id,
      -888::int AS source,
      -999::int AS target,
      cost::float * ' || source.prec || ' * ' || target.prec || ' AS cost
    FROM n50.n50_vegsti
    WHERE ogc_fid = ' || source.id || '

    UNION SELECT
      ogc_fid AS id,
      -999::int AS source,
      -888::int AS target,
      cost::float * ' || target.prec || ' * ' || source.prec || ' AS cost
    FROM n50.n50_vegsti
    WHERE ogc_fid = ' || target.id;
  END IF;

  -- This is the actual SQL query that takes the edges and applies the
  -- K-Shortest Path routing algorithm implemeted by `pgr_ksp`.

  sql := 'SELECT
    path.path_id,
    ST_LineMerge(ST_Collect(vegsti.geometri)) AS geom,
    SUM(path.cost) AS cost
  FROM
    pgr_ksp(
      ''' || sql || ''',
      -888,
      -999,
      ' || targets || ',
      directed:=false
    ) AS path,
    n50.n50_vegsti AS vegsti
  WHERE path.edge = vegsti.ogc_fid
  GROUP BY path.path_id
  ORDER BY SUM(path.cost)';

  -- Now, get the shortest paths and process them before returning

  FOR rec IN EXECUTE sql
  LOOP
    -- Return if no route was found
    IF rec.geom IS NULL THEN
      RAISE NOTICE '[ROUTER] route geometry is NULL; returning';
      RETURN;
    END IF;

    -- If the geometry is not a `LINESTRING` it means `ST_LineMerge` was not
    -- able to merge path geometry to a contiguous line.
    IF GeometryType(rec.geom) != 'LINESTRING' THEN
      RAISE NOTICE '[ROUTER] route geometry is MULTILINE; returning';
      RETURN;
    END IF;

    -- Ok, so `rec.geom` is the combined geometry of the all edges in the
    -- shortes path. However, since we rarly route from the exact beginning of
    -- an edge we need to trim the route in both ends to the points where we
    -- want the route to source and target.

    -- Locate the point on the route closest to the source point
    prec1 := ST_LineLocatePoint(rec.geom, point1);

    -- Locate the point on the route closest to the target point
    prec2 := ST_LineLocatePoint(rec.geom, point2);

    RAISE NOTICE '[ROUTING] path_id=%, start=%, end=%', rec.path_id, prec1, prec2;

    -- Now that we have the closest points in both ends we just need to trim the
    -- route in order to fit our source and target requirements. Since the route
    -- can be returned source to target or target to source we check reverse the
    -- route when necessary before trimming.

    IF (prec2 < prec1) THEN
      rec.geom := ST_Reverse(ST_LineSubstring(rec.geom, prec2, prec1));
    ELSE
      rec.geom := ST_LineSubstring(rec.geom, prec1, prec2);
    END IF;

    RAISE NOTICE '[ROUTING] path_id=%, cost=%', rec.path_id, rec.cost;

    -- Return the resulting record as defined in our function definition
    path_id := rec.path_id;
    cost := rec.cost;
    geom := rec.geom;

    -- Continue processing the next shortest route
    RETURN NEXT;

  END LOOP;
  RETURN;

END;

$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;