Selecting touching polygons by PostGIS for QGIS atlas?

Question

I‘ve a PostGIS polygon layer where each feature represents a page of an atlas in the QGIS print composer.

CREATE TABLE public.composer (
gid serial NOT NULL,
geom geometry(Polygon, 31256),
page_number varchar(8),
page_north varchar(8),
page_west varchar(8),
page_south varchar(8),
page_east varchar(8),
CONSTRAINT composer_pkey PRIMARY KEY (gid));
;

INSERT INTO public.composer (geom, page_number) VALUES
(ST_GeometryFromText('Polygon((0 0, 300 0, 300 200, 0 200, 0 0))', 31256), 'A-A'),
(ST_GeometryFromText('Polygon((300 0, 600 0, 600 200, 300 200, 300 0))', 31256), 'A-B'),
(ST_GeometryFromText('Polygon((0 0, 300 0, 300 -200, 0 -200, 0 0))', 31256), 'A-C'),
(ST_GeometryFromText('Polygon((300 0, 600 0, 600 -200, 300 -200, 300 0))', 31256), 'A-D')
;

What I want to do is insert the numbers of the touching features into the attributes of each feature.

Can anyone help me with the SQL query?

Doing this manually is quite boring since my atlas contains around 100 features.

duplicate of https://gis.stackexchange.com/questions/214300/how-to-determine-neighbouring-tile-ids-in-qgis ? — Ian Turton, Dec 05 '16 at 08:59
I will try Carrillo's processing script within the next few days. However, it would be nice to find a PostGIS query for this. — eclipsed_by_the_moon, Dec 11 '16 at 22:25
So could you tell us what your definition of touches is? What do you want inserted for each of those rows you've provided? — Evan Carroll, Dec 12 '16 at 03:27

score 3 · Accepted Answer · answered Dec 12 '16 at 00:21

This is a bit tricky, because your polygons aren't normalized in terms of orientation or starting point: some are oriented counter-clockwise, starting in the lower-left, while other are oriented clockwise, starting in the upper-left. This means that you can't use simple index-based calls to ST_PointN to pull out the the edge representing a specific edge of a polygon.

How to fix this? I'd run a simple

UPDATE composer SET geom = ST_Envelope(geom);

Now all of your polygons will be oriented clockwise, starting in the lower-left. With these standardized polygons, we know that:

the western edge is defined by points 0 and 1
the northern edge is defined by points 1 and 2
the eastern edge is defined by points 2 and 3
the southern edge is defined by points 3 and 4

Using this, you can use ST_PointN to extract the two endpoints that define each of the edges. A couple "gotchas": ST_PointN uses indexes that start at 1, not 0, and it only runs on LineStrings, not polygons. So, we use ST_Boundary to get the line representing our polygon's exterior ring. Now, we can pull out the endpoints for each edge and join the table back to itself:

UPDATE composer c
SET page_west = 
    (SELECT page_number FROM composer west 
     WHERE west.gid != c.gid
     AND   ST_Intersects(west.geom, ST_PointN(ST_Boundary(c.geom), 1))
     AND   ST_Intersects(west.geom, ST_PointN(ST_Boundary(c.geom), 2))),
    page_north =  
    (SELECT page_number FROM composer north
     WHERE north.gid != c.gid 
     AND   ST_Intersects(north.geom, ST_PointN(ST_Boundary(c.geom), 2))
     AND   ST_Intersects(north.geom, ST_PointN(ST_Boundary(c.geom), 3))),
    page_east =  
    (SELECT page_number FROM composer east 
     WHERE east.gid != c.gid 
     AND   ST_Intersects(east.geom, ST_PointN(ST_Boundary(c.geom), 3))
     AND   ST_Intersects(east.geom, ST_PointN(ST_Boundary(c.geom), 4))),
    page_south =  
    (SELECT page_number FROM composer south 
     WHERE south.gid != c.gid         
     AND   ST_Intersects(south.geom, ST_PointN(ST_Boundary(c.geom), 4))
     AND   ST_Intersects(south.geom, ST_PointN(ST_Boundary(c.geom), 5)));

Selecting touching polygons by PostGIS for QGIS atlas?

1 Answers1