[postgis-users] Dot Density idea

[John Abraham]

Just to follow up on this, here is the code we (Abdel-Rahman M. Muhsen and I) finally got around to writing.  We found a function called randompoint which just creates random points within the bounding box of the geometry from Alexandre Sorokine. http://sorokine.blogspot.com/2011/05/postgis-function-for-random-point.html (I see he's updated it 4 days ago to be better for multi-polygons!  We don't have his updates yet.)

Then we wrote a function called dot_density which creates a table called dp which has a point geometry for each point in the polygon.  If there are negative numbers in the point count, it creates the positive number of points but flags them in another column. (We often plot dot density maps with red dots for decreases and blue dots for increases.)

There are still a few improvements we'll likely make.  We'd like to create the geometry column using the addgeometrycolumn function.  We'd like to specify the output table name as a function parameter.  There was also a good idea to use Halton sequences or other "pseudo-random" sequences, instead of truly random points, and Martin Davis implemented a few of the other "pseudo-random" ideas in JTS http://lin-ear-th-inking.blogspot.com/2010/05/more-random-points-in-jts.html .  Other suggestions are welcomed.

This is how we typically use it:

create or replace view dot_cnt as 
select the_geom, polygon_id, (quantity_column/1000000)::integer as numpoints from original_data;

select dot_density('dot_cnt',
'the_geom',
'polygon_id',
'numpoints');

alter table dp rename to save_it_for_future_use;

And here are the functions:


CREATE OR REPLACE FUNCTION randompoint(geom geometry)
  RETURNS geometry AS
-- from Alexandre Sorokine
$BODY$
DECLARE
 maxiter INTEGER := 1000;
 i INTEGER := 0;
 x0 DOUBLE PRECISION;
 dx DOUBLE PRECISION;
 y0 DOUBLE PRECISION;
 dy DOUBLE PRECISION;
 xp DOUBLE PRECISION;
 yp DOUBLE PRECISION;
 rpoint Geometry;
BEGIN
 -- find envelope
 x0 = ST_XMin(geom);
 dx = (ST_XMax(geom) - x0);
 y0 = ST_YMin(geom);
 dy = (ST_YMax(geom) - y0);
 
 WHILE i < maxiter LOOP
  i = i + 1;
  xp = x0 + dx * random();
  yp = y0 + dy * random();
  rpoint = ST_SetSRID( ST_MakePoint( xp, yp ), ST_SRID(geom) );
  EXIT WHEN ST_Within( rpoint, geom );
 END LOOP;
 
 IF i > maxiter THEN
  RAISE NOTICE 'number of interations exceeded max';
 END IF; 
 
 RETURN rpoint;
END; 
$BODY$
  LANGUAGE 'plpgsql' VOLATILE
  COST 100;
ALTER FUNCTION randompoint(geometry) OWNER TO "usrPostgres";

-- Function: dot_density(text, text, text, text)

CREATE OR REPLACE FUNCTION dot_density(geom_table text, geom_col text, zone_col text, num_of_points_col text)
  RETURNS SETOF record AS
$BODY$
DECLARE 
 counter integer:=0;
 tazrec record;
 pointrec record;
 result record;

 num_points integer:=0;
 np integer :=0; 
begin

DROP sequence if exists randpnt_id;
CREATE sequence randpnt_id;

DROP TABLE IF EXISTS dp;
CREATE TABLE dp
(
  gid integer PRIMARY KEY,
  ser integer,
  "zone" integer,
  decrease_or_increase integer,
  the_geom geometry
); 

for tazrec in EXECUTE 'SELECT ' || zone_col || ' as geom_col , ' || zone_col || ' as zone_col, '|| num_of_points_col || ' as num_of_points_col FROM ' || geom_table Loop
 RAISE INFO 'Treating zone: %' , tazrec.zone_col;
 num_points = tazrec.num_of_points_col;

 IF num_points !=0 THEN np := num_points/abs(num_points);
 ELSE np=0;
 END IF;
 
 

 EXECUTE 'INSERT INTO dp SELECT nextval(''randpnt_id'') as gid, generate_series, '|| tazrec.zone_col ||  ', ' || np ||' , randompoint(the_geom) FROM ' || geom_table || ', generate_series(1, '|| abs(num_points) ||') WHERE '|| zone_col || '='|| tazrec.zone_col ;
  

END LOOP;

RETURN ; 
end;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;
ALTER FUNCTION dot_density(text, text, text, text) OWNER TO postgres;


On 2010-05-06, at 3:10 AM, strk wrote:

> ST_RandomPoinsOnSurface(geometry, numpoints) would be an interesting
> function indeed. Sounds like a good job for GEOS/JTS.
> 
> --strk;
> 
> On Mon, May 03, 2010 at 10:49:32PM -0600, John Abraham wrote:
>> One of the things I miss about using ESRI's GIS is the ability to do dot-density maps.  Within a polygon, the number of dots is proportional to a value, and the dots are randomly placed.  I find it useful to be able to present several data values at once (e.g. blue dots for population, red dots for employment).  
>> 
>> I also find that it is a more intuitive way of scaling for zone size than dividing the value by the area of the zone.  That is, the count of the dots represents the actual number, but the density of the dots represents the density of the number.  So I don't have to decide whether to divide the value by the area of the polygon to plot density: both the absolute number and the density are easily visible.
>> 
>> Since my open-source GIS viewing systems (mostly QGIS and Mapserver) won't plot dot-density, I've done without.
>> 
>> But today I realized that I can build these on the server instead.  I can generate random points within the bounding-box of the polygon, throwing out those that aren't contained within the polygon, repeating until I have enough.  Then I can save these points as a separate layer, and display this layer using almost any desktop or web based viewer!
>> 
>> Has anyone done this?  Can I do it in SQL or do I need to write something in PL/pgsql?
>> 
>> --
>> John Abraham
>> 
>> PS I just bought the Postgis In Action book; enjoying it so far.
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 
> -- 
> 
>  ()   Free GIS & Flash consultant/developer
>  /\   http://strk.keybit.net/services.html
> 
> 

-- 
John Abraham
This email address is being phased out, please update your records to use john at theabrahams.ca for personal email and jea at hbaspecto.com for work email.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20110527/36d7fd02/attachment.html>