[GIS] Determine which points are contained inside a given polygon

point-in-polygonpolygonpostgispostgresqlsql

I have these database tables defined:

CREATE TABLE regions (
  id SERIAL PRIMARY KEY,
  name VARCHAR(150) NOT NULL
)

CREATE TABLE region_points (
  id SERIAL PRIMARY KEY,
  region_id SERIAL REFERENCES regions (id) ON DELETE CASCADE,
  index INTEGER NOT NULL,
  point GEOMETRY(Point, 4326) NOT NULL
)

CREATE TABLE stations (
  id SERIAL PRIMARY KEY,
  name VARCHAR(150) NOT NULL,
  point GEOMETRY(Point, 4326) NOT NULL
)

A region is a polygon, its points being defined in the region_points table. A station may reside at a single point somewhere inside a region polygon.

I'd like to write a single query that gets all stations inside a particular region polygon. At the moment, I don't know how, so I do it in multiple steps:

# 1. Get all points that comprise the region polygon.
SELECT point
FROM region_points
WHERE region_id = 1234
ORDER BY index

# 2. Programmatically (non-SQL) construct a linestring from the points. Not shown here.

# 3. Get all stations inside the polygon.
SELECT id, name, point
FROM (
    SELECT ST_GeomFromText($1, 4326) AS polygon
) AS query1, (
    SELECT id, name, point
    FROM stations
) AS query2
WHERE ST_Contains(query1.polygon, query2.point)

As you can see, I have to construct the linestring programmatically and then pass it into the query as an argument to ST_GeomFromText, because I don't know how to do it in SQL.

I don't want to do anything programmatically. How can I re-write all of this as one SQL query?

Best Answer

The general workflow to construct your polygon would be something like this: st_makepolygon(st_linefromMultipoint(st_collect(points.geom)))

But you have to be aware that :

  • your linestring has to be closed
  • order of vertices does mater

a rudimary sql querry to create a polygon from a set of points would be the following: 

with line as (
with mp as (
SELECT 
-- multipoint
ST_Collect(ARRAY[
ST_GeomFromText('POINT(1 2)'),
ST_GeomFromText('POINT(2 3)'),
ST_GeomFromText('POINT(-2 3.2)'),
ST_GeomFromText('POINT(0 0)')
]) as geom -- from point table where = [exrp] order by
)
select (ST_AddPoint( st_linefromMultipoint(mp.geom),ST_PointN(st_linefromMultipoint(mp.geom),1) )) as geom from mp
)
select st_makepolygon(line.geom) from line

Afterwards since you have constructed your polygon just wrap the above in your main query: 

with polygon as (
-- above query 
)
select count(a.fid) from stations a, polygon b where st_contains(b.geom,a.point)

ofc it is very possible that your nodes in the region table do not follow an ordering at all you could use ST_ConvexHull functions to 'simplify' the polygon creation:

with mp as (
SELECT 
-- multipoint mu
ST_Collect(ARRAY[
ST_GeomFromText('POINT(1 2)'),
ST_GeomFromText('POINT(2 3)'),
ST_GeomFromText('POINT(-2 3.2)'),
ST_GeomFromText('POINT(0 0)')
]) as geom -- from point table where = [exrp] order by
)
select st_convexhull(mp.geom) as geom from mp
Related Question