5

I have a number of point and line tables in postgres in which I want to update an attribute according to which area of a 5km grid they fall into. I use the following postgis query to do this:

UPDATE point_table AS B
SET tile = A.tile_name
FROM 5km_grid_polygons AS A
WHERE ST_Within (B.geom, A.geom);

This worked find for most of my data. However, when a point or line falls exactly on a grid boundary (as in screenshot below), then no attribute is copied over.

enter image description here

This has happened hundreds of times, so manually fixing is not an option. I am wondering if anyone can think of a way to create a rule that assigns one grid square to the point or line that is 'in limbo'?

The rule/query would have to be consistent, for instance, always assigning the attribute from the lower or left grid square to the point or line. It cannot pick one of the adjacent grid squares at random.

4

For each point, you can find all intersecting grid cells using a subquery, and then sort them so that you get a consistent assignment. The example below would give preference to the lower-left grid cell:

UPDATE point_table
SET tile = (SELECT tile_name 
            FROM 5km_grid_polygons
            WHERE ST_Intersects(point_table.geom, 5km_grid_polygons.geom)
            ORDER BY ST_XMin(5km_grid_polygons.geom),
                     ST_YMin(5km_grid_polygons.geom)
            LIMIT 1);
| improve this answer | |
2

The semantics of "within" require the interior of the polygon to intersect the point, but the problem points are exactly on the boundary. See DE-9IM for the formal definition of these things. A better spatial predicate to find intersections on both exterior and interior is "intersects" (i.e. ST_Intersects).

A point on a shared border needs to be matched to only one grid record, which can be done with a SELECT DISTINCT clause. Choosing a bias towards picking a grid name from the upper-right side would look like this:

UPDATE point_table AS A
SET tile = B.tile_name
FROM (
  SELECT DISTINCT ON (p.id)
    p.id, g.tile_name
  FROM point_table p
  JOIN "5km_grid_polygons" g ON ST_Intersects(p.geom, g.geom)
  ORDER BY p.id, ST_XMin(g.geom) DESC, ST_YMin(g.geom) DESC
) B
WHERE B.id=A.id;

where it is assumed you have a primary key id on the point table to make the subquery work.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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