3

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

enter image description here

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.

3
  • 1
    duplicate of gis.stackexchange.com/questions/214300/… ?
    – Ian Turton
    Dec 5, 2016 at 8: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. Dec 11, 2016 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? Dec 12, 2016 at 3:27

1 Answer 1

3
+100

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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