Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Hi all hoping someone can help me out with this issue.

I've got a table of polygons which are grid squares which get overlayed on a map. These are used when querying data so that points get grouped into each grid square. All this works fine till a point is on the boundary between the grid squares. I've included an example of this below. If I use ST_Contains I get no results as the point isn't contained within any of the polygons. If I use ST_Intersects the point intersects with the boundaries of multiple grid squares which isn't suitable.

Is there a way I can ensure that the point is only considered to be in one polygon? In are previous system before migrating to using postgis we used the PNPOLY algorithm that ensured this.

DROP TABLE IF EXISTS areasTemp;
    CREATE TEMP TABLE areasTemp
    (
      id integer,
      areacode character varying(10),
      polygon geometry
    );
    INSERT INTO areasTemp VALUES (634,'1x50',ST_GeomFromText('POLYGON((1 50,2 50,2 51,1 51,1 50))', 4326)),
    (663,'0x51',ST_GeomFromText('POLYGON((0 51,1 51,1 52,0 52,0 51))', 4326)),
    (633,'1x51',ST_GeomFromText('POLYGON((1 51,2 51,2 52,1 52,1 51))', 4326)),
    (664,'0x50',ST_GeomFromText('POLYGON((0 50,1 50,1 51,0 51,0 50))', 4326));

-

SELECT * FROM areasTemp WHERE ST_Contains(ST_GeomFromText('POINT(1 51)', 4326), polygon);

0 rows retrieved

SELECT * FROM areasTemp WHERE ST_Intersects(ST_GeomFromText('POINT(1 51)', 4326), polygon);

4 rows retrieved

share|improve this question
    
Which single polygon were you expecting? I don't know how PNPOLY can find only one polygon from this situation. Whether it is based on some spatial preference or rule, or if it is the first record (which could be random). I'd investigate the behavour of the logic for this boundary case in your PNPOLY system. –  Mike T Jan 26 '12 at 23:49
add comment

2 Answers

I've just had the same issue when trying to count the number of fire events per 0.25 degree grid cell. Using st_intersects(fire.geom, grid.geom) in a spatial join double-counted any fires "sitting on the fence" between two grid cells.

It is however possible to distinguish which fires have been double counted. I used a PostgreSQL window function, row_number(), partitioning over the point identifier (in my case fires.id) you can just select the first row in each partition (the partition is the set of fire-gridcell intersections which are actually the same fire).

In my example below I simply select those rows with row_number=1 which chooses the point-gridcell intersection with the lowest grid cell id, due to the ordering in the partition.

select * from
(
select
st_centroid(j25.geom) as centroid,
row_number() over (partition by fires.id order by j25.id)
from fires join grid_cells
on st_intersects(grid_cells.geom, fires.geom) and 
grid_cells.geom && fires.geom
) as fire_grid_intersection
where row_number = 1
share|improve this answer
add comment

If you don't care which polygon the point is in, then you just need to add the LIMIT keyword:

SELECT * FROM areasTemp WHERE ST_Intersects(ST_GeomFromText('POINT(1 51)', 4326), polygon) LIMIT 1;

Which will choose whichever row PostgreSQL chooses first. This will be arbitrary but consistent.

If you need a more specific rule as to which polygon a point is in, then you'll need to add an ORDER BY part to the SELECT statement to ensure your preferred table row is at the top of the results list.

Edit:

Following on from your reply, you'll need to select your polygon as a subquery, and join against that:

SELECT points.*, grid.*
    FROM (SELECT * FROM grid WHERE ST_Intersects(points.the_geom) LIMIT 1) AS grid
INNER JOIN
    ... your join clause ...

Aplogies if the details are little sketchy, I'm not near a PostGIS database at the moment to play around with it, but a subquery would be the way to go.

An alternative might be to rewrite the PNPOLY function as a PL/pgSQL function and put that into your WHERE clause.

share|improve this answer
    
Thanks for the answer but unfortunately it won't work in my case. I was trying to keep the example really simple so as to not make it to confusing, but in the actual system got a data table of ~18million rows of data which the area table joins to. So can't use any limits against it. Basically for each data point need to check against the areas table and work out what polygon it belongs too. –  Mark Davidson Jan 26 '12 at 11:25
    
Just made an edit to expand on your requirements. –  MerseyViking Jan 26 '12 at 13:25
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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