1

I am trying to do a 'point outside polygon' query for GBR with ~ 4 million points. Most of these points are inside the boundary but a handful are falling way outside the boundary.

Below steps are done to pre-process point data -

  1. Add primary key GID (ALTER TABLE point_gbr ADD COLUMN gid serial PRIMARY KEY)
  2. Add GEOM (ALTER TABLE point_gbr ADD COLUMN the_geom geography(POINT,4326))
  3. Update Geom (UPDATE point_gbr SET the_geom = ST_SetSRID(ST_MakePoint(longitude,latitude),4326))
  4. Create Spatial index (CREATE INDEX pt_idx_geom ON point_gbr USING GIST(the_geom))

For Polygon data -

Uploaded the A0 shapefile using 'Postgis3 Shapefile and DBF loader' and set the SRID to 4326. While upload I also create an index on GEOM column.

I executed the below query to find the points outside the boundary only if they fall beyond 500m, but it is still running for 6+ hours. Not sure what am I doing wrong?

SELECT 
  ST_Distance(a.the_geom::geography, b.geom::geography) as distance_metre, a.*
FROM
  point_gbr  a , polygon_gbr  b
WHERE
  a.gid  NOT IN
    (SELECT a.gid 
     FROM point_gbr  a ,  polygon_gbr  b 
     WHERE ST_Within(a.the_geom::geometry, b.geom)
    )
 AND ST_Distance(a.the_geom::geography, b.geom::geography) > '500' 

The laptop has a RAM of 24 GB and 4 cores and is running PostgreSQL 11 with PostGIS 2.5.

Also, is there a better way to do this for similar complex countries like Russia, Mexico?

  • 2
    Try storing the points as geography to save the casts and subdivide the polygon and index it. – Ian Turton May 26 at 18:40
  • 2
    I believe that the first thing to do is to subdivide your geometries postgis.net/docs/ST_Subdivide.html. Then I guess I would create a new indexed column for the points "is_in" and update the value with the id of the polygon that the point is in. That way you could discard most points without going to spatial queries and for the rest you could operate with small geometries. – user30184 May 26 at 19:33
  • Thanks, I will work on st_subdivide and update on the progress. – Pankaj May 27 at 10:53
2

Several modifications are required to make it work.

Before that, let's note that geom is a confusing name for a geography field... call it geog and so it will be obvious when a cast to geography is needed or not.

A spatial index is used only if the same expression as the one used to create the index is used. In your case, it is a index on a geography field but the expression is using geom::geometry so it can't be used. Either use the geog field, or create an index on the geometry CREATE INDEX pt_idx_geom ON point_gbr USING GIST(geometry(the_geog))

Having one huge polygon is a performance killer for point in polygon operations. Create a new table with the big polygons divided in smaller ones. You can read more here.

Your points are already geography, so the simplest is to also save the subdivided polygons as geography.

Then make sure to use a function that do rely on the indexes, such as st_DWithin.

You would join the point and the polygon based on proximity, and keep only the not-joined points.

Finally, do a cross join lateral to find the nearest polygon and its distance.

It should look like:

SELECT 
  ST_Distance(a.the_geog, closest_poly.geog) as distance_metre, a.*
FROM
  point_gbr  a 
  LEFT JOIN polygon_gbr_subdivided  b
     ON ST_DWithin(a.the_geog, b.geog,500)
  CROSS JOIN LATERAL
   (SELECT
     a.geom <-> b2.geom as dist
     FROM polygon_gbr_subdivided  b2
     ORDER BY a.geog <-> b2.geog
   LIMIT 1) AS closest_poly
WHERE b.id IS NULL;

Now this query may be tweaked should the lateral join be executed before the left join filtering... only a explain (buffers,analyze) ... will tell you.

| improve this answer | |
  • Will there be issues with the transformed subdivided polygons not "lining up" along the shared boundaries? – dr_jts May 26 at 20:01
  • @dr_jts I don't think so because there is no reprojection. The subdivided shapes should be topologically correct – JGH May 26 at 20:13
  • Thanks JGH, that was helpful. I will be implementing the suggested changes and will post on the improvement. – Pankaj May 27 at 10:53

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.