Points in polygon by months

Question

I have two tables: points(gid,id,datum,x,y,geom) and hexa_grid(gid,geom).

I must create new table hexa_grid_pip (gid,2008_01,2008_02 ... 2011_12, geom) In column 2008_01 ... 2011_12 must be number points in polygon. Columns 2008_01 ... 2011_12 are shortcuts of months.

Is it possible to do this by using script in PostGIS, or do I have to count month by month in QGIS?

Answer 1

Here's a partial answer. To begin with, Postgresql cannot have column names beginning with a digit. So you'll have to name the year_month columns differently. It could be "_2008_01" (note the leading underscore) or something like ym_2008_01. I'll assume the second naming convention.

First create the hexa_grid_pip table:

CREATE TABLE hexa_grid_pip AS 
  SELECT gid FROM hexa_grid;

This will make the new table with only one column, gid, using the gid values from the hexa_grid table. Now add all the other year_month columns:

ALTER TABLE hexa_grid_pip ADD COLUMN ym_2008_01 INTEGER, 
  ADD COLUMN ym_2008_02 INTEGER, 
  ADD COLUMN ym_2008_03 INTEGER, 
  ADD COLUMN ym_2008_04 INTEGER,
  ....
  ;

Now, you can update one of the month columns with an expression like:

UPDATE hexa_grid_pip SET ym_2008_01=( 
  SELECT COUNT(p.gid) 
  FROM points AS p, hexa_grid AS h
  WHERE ST_Within(p.geom, h.geom) AND
    date_part('year', p.datum) = 2008 AND
    data_part('month',p.datum) = 1 AND
    h.gid = hexa_grid_pip.gid
  );

Now, what you'll want to do is wrap the above into a function which does two nested loops, over the 4 years and 12 months (this is just a skeleton example, not tested!):

DECLARE
yr integer :=2008;
mo integer :=1;
col text;

BEGIN
FOR yr IN 2008..2011 LOOP
  FOR mo IN 1..12 LOOP

  -- Make the column header
  col='ym_' || to_char(yr, '9999') || '_' || to_char(mo,'99');

  -- do the UPDATE similar to above, but with the yr, mo and col variables
      UPDATE hexa_grid_pip SET col=( 
      SELECT COUNT(p.gid) 
      FROM points AS p, hexa_grid AS h
      WHERE ST_Within(p.geom, h.geom) AND
        date_part('year', p.datum) = yr AND
        data_part('month',p.datum) = mo AND
        h.gid = hexa_grid_pip.gid;

  END LOOP;
END LOOP;
END;

Again, I've not tested the above. It will require debugging. HTH

Micha

Answer 2

I think you should divide this in two separate things.

First you have the calculations to find how many points you have in every polygon each month.

The other part is the presentation, or report part, where you make a pivot table out of it. This is because, what you are doing is creating a huge table, that will grow in width every month.

From a database perspective is better to handle the result vertical like: hexa_grid_pip(gid, year_month, number_of_points)

Then, when someone wants to look at the data you create a pivot table of the months that is interesting, or visualize the result in any other way. The point is that that part is about presentation, not how to get the data.

To get the data you can do something like:

SELECT count(*) AS number_of_points,the_year, the_month, hid
FROM
(
    SELECT 
        EXTRACT(month from datum) as the_month, 
        EXTRACT(year from '2012-01-01'::timestamp) AS the_year,
        p.geom AS the_point, 
        h.geom AS the_polygon, 
        p.gid AS pid, 
        h.gid AS hid
    FROM 
        points p INNER JOIN hexa_grid h 
        ON ST_Intersects(p.geom, h.geom)
) AS pip
GROUP BY hid, the_year, the_month;

Then to make a pivot you can use function "crosstab" or something like that:http://www.postgresql.org/docs/9.1/static/tablefunc.html