2

I am a new user in a PostGIS environment and don't even know where to start.

I have two tables, one storing lines and the other storing polygons. I to create a query that can be used to get a separate table where each line has a separate field listing the IDs of the polygons that the line touches. The polygons are represented as a grid of squares. And the lines are a grid of streets. Green squares are polygons, others are lines

As a result, we would like to get the following list

result table

1
  • You are probably looking for a cross lateral join. Can you provide some example data? I will be able to provide an answer then. May 16 at 9:00

1 Answer 1

4

This is generally solved using an [INNER] JOIN, or if you are also interested in lines that do not intersect any grid cells in the result set, a LEFT JOIN; you then want to GROUP BY your line_ids and collect the name_polygon with an ARRAY_AGG - the spatial predicate function you are looking for is ST_Intersects:

SELECT
  l.line_id,
  ARRAY_AGG(DISTINCT p.polygon_name) AS polygon_names
FROM
  <lines> AS l
  JOIN
  <polygon> AS p
    ON ST_Intersects(l.geom, p.geom)
GROUP BY
  l.line_id
;

to only return lines that actually intersect at least one polygon, or

SELECT
  l.line_id,
  ARRAY_AGG(DISTINCT p.polygon_name) AS polygon_names
FROM
  <lines> AS l
  LEFT JOIN
  <polygon> AS p
    ON ST_Intersects(l.geom, p.geom)
GROUP BY
  l.line_id
;

to also include those that do not intersect any polygons (receiving NULL as polygon_names).


Considerations:

  • if you expect that the majority (or all) rows in <lines> will get returned, and especially

    1. if you need to return more than one column, or others than those covered by the PRIMARY KEY or any multi-column index
    2. if you need to apply transformations to any of the returned values from <polygon>

    it is more performant to run a [LEFT|CROSS JOIN] LATERAL to avoid the costly extended GROUP BY list:

    SELECT
      l.line_id,
      g.polygon_names
    FROM
      <lines> AS l
      CROSS JOIN LATERAL (
        SELECT
          ARRAY_AGG(DISTINCT p.polygon_name) AS polygon_names
        FROM
          <polygon> AS p
        WHERE
          ST_Intersects(l.geom, p.geom)
      ) AS g
    ;
    

    to only return lines that actually intersect at least one polygon, or

    SELECT
      l.line_id,
      g.polygon_names
    FROM
      <lines> AS l
      LEFT JOIN LATERAL (
        SELECT
          ARRAY_AGG(DISTINCT p.polygon_name) AS polygon_names
        FROM
          <polygon> AS p
        WHERE
          ST_Intersects(l.geom, p.geom)
      ) AS g ON TRUE
    ;
    

    to also include those that do not intersect any polygons (receiving NULL as polygon_names)

  • ST_Intersects (as well as most other 2D spatial predicate functions) are running against the GEOS backend, which has recently been enhanced with prepared geometry caching that boosts the underlying computations; this has a massive effect if you expect to have many overlaps with a (few) potentially large and complex geometries each.

    In all others cases it is worth comparing performance to the native and highly optimized ST_DWithin function with 0 distance - simply replace ST_Intersects(l.geom, p.geom) with ST_DWithin(l.geom, p.geom, 0) in the above queries - on a test with 100k random segments over 1M grid cells, the factor of execution time is almost 3:1 in favor of ST_DWithin


Run any of the above queries within a CREATE TABLE statement:

CREATE TABLE <name> AS (
  <query>
);

to get a new table with the results.

1
  • 2
    Note that perfromance depends greatly on the presence of a spatial (e.g. GIST) index on the <polygon> geometry column!
    – geozelot
    May 16 at 11:13

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