Points in polygon by months

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?

• Where does the month data come from? I don't see any dates in your points table? – Micha Jul 25 '12 at 12:50
• At least in Swedish, (and other nordic languages) "datum" means date. – Nicklas Avén Jul 26 '12 at 8:22
• In Czech languages "datum" is "date" in English. Column "datum" is type timestamp (2008-01-09 16:38:00.000) – user9128 Jul 26 '12 at 9:13
• Ah, OK. (In GIS, datum refers to something else...). With this datum column as date, I put together a rough answer below. – Micha Jul 26 '12 at 11:59

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,
....
;


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

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

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