Determine which points are contained inside a given polygon

Question

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?

Answer 1

I'd change my tables to something like this:

CREATE TABLE regions (
  id SERIAL PRIMARY KEY,
  name text NOT NULL,
  geom geometry(Polygon,4326)
)

CREATE TABLE stations (
  id SERIAL PRIMARY KEY,
  region_id INTEGER NOT NULL REFERENCES regions (id) ON DELETE CASCADE,
  name text NOT NULL,
  point GEOMETRY(Point, 4326) NOT NULL,
  CONSTRAINT fk_reg FOREIGN KEY (region_id)
    REFERENCES regions (id) MATCH SIMPLE
    ON UPDATE CASCADE ON DELETE CASCADE
)

Notes:

• you can use a polygon geometry
• text is prefered to varchar in 99% in PostgreSQL
• a foreign key is better than doing any spatial queries (about 1000x faster, more error-proof and less work)

If you just need to select stations in a region:

SELECT id, name, point
FROM stations 
WHERE region_id = 123;

Keep this structure if your goal is to have a useful database. Keeping points and making polygons in a query is like hunting ants with a tank - it's foolish unless you need the points much more often than the polygons (the questions suggest otherwise, but we don't know how much did you simplify the problem).

EDIT: If your don't know which station lies in which region, you can use some other tool (and save overhead with creating the database etc.), or you can compute it using PostGIS.

If you want to use PostGIS for this, then first step is to create and populate the regions table. If you have coordinates, make the polygons using WKT; if you have points somewhere in PostGIS, call ST_MakePolygon. Or use Nickves' answer for this. Don't create extra table for points unless you already have it and too much other features depend on it. If the points are used about as much as the polygon, ST_DumpPoints is your friend; extracting points should be faster than making a polygon.

The second step would be to prepare the stations table. The only difference from above is that you omit the NOT NULL for the first time (you can add it faster to make it clean). Create an AFTER INSERT trigger on Stations table calling following function:

CREATE OR REPLACE FUNCTION create_region_id() RETURNS trigger AS $$
BEGIN
  IF NEW.region_id IS NULL THEN
    UPDATE stations
    SET regions_id = (
      SELECT id FROM region 
      WHERE ST_CONTAINS(region.geom, NEW.point)
    )
    WHERE id = NEW.id;
  END IF;
  RETURN NEW;
END;
$$ LANGUAGE plpgsql;

Then populate your stations table whatever way you want, leaving region_id NULL.

EDIT2: for the first time I misread your queries and thought you oversimplified things much more than you did. Sorry.

Answer 2

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