[postgis-users] fishnet for counting points
Stephen Woodbridge
woodbri at swoodbridge.com
Mon Oct 17 17:02:39 PDT 2011
Matthias,
This is a different approach but it might be useful. This is untested,
but hopefully it will give you the idea:
select
floor((st_x(a.the_geom)-b.minx)/1000)::integer as col,
floor((st_y(a.the_geom)-b.miny)/1000)::integer as row,
count(*) as cnt
from
(
select floor(min(st_x(the_geom)))::integer as minx,
floor(min(st_y(the_geom)))::integer as miny
from points ) as b,
points a
group by col, row
order by cnt desc;
With the col, row, values and the count you can create the polygons
after the fact.
-Steve W
On 10/17/2011 6:07 PM, Matthias Ludwig wrote:
> Hello,
>
> for a statistical analysis I need to create a polygon grid to count the points inside. The analysis area measures about 15x15km. At the moment I use the code below. It's functioning but is slow for high point counts(>5000-100000 points) and small polygon sizes (5-20m). Because the points are clustered and the big area size it's not possible to insert all polygons and delete those which are empty. Presently one polygon is inserted for every point in the first step. In the second step dublicated polygons where deleted. How can I insert one polygon in a specific place?
>
> CREATE TABLE point_fishnet (gid serial not null primary key, area text, count integer default 0);
> SELECT addgeometrycolumn('point_fishnet','the_geom', 900913, 'POLYGON', 2);
>
> INSERT INTO point_fishnet(the_geom, area)
> (
> SELECT ST_Translate(bb.the_geom, x_koord, y_koord) AS the_geom,
> 'germany'
> FROM generate_series(
> (
> SELECT floor(min(ST_X(the_geom)))::integer
> FROM points
> ),
> (
> SELECT ceil(max(ST_X(the_geom)))::integer
> FROM points
> ),
> 1000 --Offset = Polygonsize
> ) AS x_koord,
> generate_series(
> (
> SELECT floor(min(ST_Y(the_geom)))::integer
> FROM points
> ),
> (
> SELECT ceil(max(ST_Y(the_geom)))::integer
> FROM points
> ),
> 1000 --Offset = Polygonsize
> ) AS y_koord,
> (
> SELECT ST_MakeEnvelope(0, 0, 1000, 1000, 900913) AS the_geom
> ) AS bb,
> (
> SELECT the_geom
> FROM points
> ) AS point
> WHERE ST_Intersects(
> ST_Translate(bb.the_geom, x_koord, y_koord),
> point.the_geom
> )
> );
>
> -- delete dublicates
> DELETE FROM point_fishnet
> WHERE gid NOT IN
> (
> SELECT min(gid)
> FROM point_fishnet
> GROUP BY the_geom HAVING count(*)>= 1
> );
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
More information about the postgis-users
mailing list