From b8fe05b388d27b20f50ec76b6fe88feb8dd2e745 Mon Sep 17 00:00:00 2001 From: Stuart Lynn Date: Mon, 3 Oct 2016 21:22:40 +0000 Subject: [PATCH] CDB_DasymetricDotDensity and fixed for basic DD --- src/pg/sql/30_dot_density.sql | 203 +++++++++++++++++++++++++++------- 1 file changed, 163 insertions(+), 40 deletions(-) diff --git a/src/pg/sql/30_dot_density.sql b/src/pg/sql/30_dot_density.sql index b0f3604..dae4ff1 100644 --- a/src/pg/sql/30_dot_density.sql +++ b/src/pg/sql/30_dot_density.sql @@ -10,45 +10,168 @@ -- misses per point the funciton accepts before giving up. -- -- Returns: Multipoint with the requested points -CREATE OR REPLACE FUNCTION cdb_dot_density(geom geometry , no_points Integer, max_iter_per_point Integer DEFAULT 1000) -RETURNS GEOMETRY AS $$ -DECLARE - extent GEOMETRY; - test_point Geometry; - width NUMERIC; - height NUMERIC; - x0 NUMERIC; - y0 NUMERIC; - xp NUMERIC; - yp NUMERIC; - no_left INTEGER; - remaining_iterations INTEGER; - points GEOMETRY[]; - bbox_line GEOMETRY; - intersection_line GEOMETRY; -BEGIN - extent := ST_Envelope(geom); - width := ST_XMax(extent) - ST_XMIN(extent); - height := ST_YMax(extent) - ST_YMIN(extent); - x0 := ST_XMin(extent); - y0 := ST_YMin(extent); - no_left := no_points; - LOOP - if(no_left=0) THEN - EXIT; - END IF; - yp = y0 + height*random(); - bbox_line = ST_MakeLine( - ST_SetSRID(ST_MakePoint(yp, x0),4326), - ST_SetSRID(ST_MakePoint(yp, x0+width),4326) - ); - intersection_line = ST_Intersection(bbox_line,geom); - test_point = ST_LineInterpolatePoint(st_makeline(st_linemerge(intersection_line)),random()); - points := points || test_point; - no_left = no_left - 1 ; - END LOOP; - RETURN ST_Collect(points); -END; +CREATE OR REPLACE FUNCTION CDB_DotDensity(g geometry, no_points integer, max_iter integer DEFAULT 1000) + RETURNS SETOF geometry +AS $$ + DECLARE + extent GEOMETRY; + test_point Geometry; + width NUMERIC; + height NUMERIC; + x0 NUMERIC; + y0 NUMERIC; + no_left INTEGER; + sample_points GEOMETRY[]; + points GEOMETRY[]; + BEGIN + + extent := ST_Envelope(g); + max_iter := 0; + width := ST_XMax(extent) - ST_XMIN(extent); + height := ST_YMax(extent) - ST_YMIN(extent); + x0 := ST_XMin(extent); + y0 := ST_YMin(extent); + no_left := no_points; + + LOOP + IF(no_left<=0 or max_iter=1000) THEN + RETURN; + END IF; + + with random_points as( + SELECT CDB_LATLNG(y0 + height*random(), x0 + width*random()) as p + FROM generate_series(1,no_left) + ) + SELECT array_agg(p) from random_points + WHERE ST_WITHIN(p, g) + into sample_points; + + RETURN QUERY select * from unnest(sample_points); + + IF sample_points IS NOT null THEN + no_left = no_left - array_length(sample_points,1); + END IF; + max_iter = max_iter + 1; + END LOOP; + + RETURN; + END +$$ LANGUAGE plpgsql; + +-- +-- Creates N points randomly distributed in the specified secondary polygons +-- +-- @param g - array of the geometries to be turned in to points +-- +-- @param no_points - the number of points to generate +-- +-- @params max_iter_per_point - the function generates points in the polygon's bounding box +-- and discards points which don't lie in the polygon. max_iter_per_point specifies how many +-- misses per point the funciton accepts before giving up. +-- +-- Returns: Multipoint with the requested points + + + +-- +-- Generate a random response based on the weights given +-- +-- @param array_ids an array of ids representing the category to return +-- +-- @param weights an array of weights for each category +-- +-- Returns : The randomly selected ID. + +CREATE OR REPLACE function _cdb_SelectRandomWeights(array_ids numeric[], weights numeric[]) returns NUMERIC + as $$ + DECLARE + result NUMERIC; + BEGIN + + WITH idw as ( + select unnest(array_ids) as id, unnest(weights) as percent + ), + CTE AS ( + SELECT random() * (SELECT SUM(percent) FROM idw) R + ) + SELECT * + FROM ( + SELECT id, SUM(percent) OVER (ORDER BY id) S, R + FROM idw as percent CROSS JOIN CTE + ) Q + WHERE S >= R + ORDER BY id + LIMIT 1 + into result; + return result; + END + $$ LANGUAGE plpgsql; + +-- +-- Weighted Dot Density +-- +-- @param no_points the number of points to generate +-- +-- @param geoms the target geometries to place the points in +-- +-- @param weights the weight for each of the target polygons +-- +-- RETURNS set of points + +CREATE OR REPLACE FUNCTION _cdb_WeightedDD(no_points numeric, geoms geometry[], weights numeric[]) + RETURNS SETOF geometry +AS $$ +DECLARE + i NUMERIC; + ids NUMERIC[]; + perGeom NUMERIC[]; + selected_poly NUMERIC; +BEGIN + with idseries as ( + select generate_series(1,array_upper(geoms,1)) as id + ) + select array_agg(id) from idseries into ids; + + FOR i in 1..no_points + LOOP + select cdb_crankshaft._cdb_SelectRandomWeights(ids, weights) INTO selected_poly; + perGeom[selected_poly] = coalesce(perGeom[selected_poly] + 1, 0 ); + END LOOP; + + raise notice 'pergeom %', perGeom; + + FOR i in 1..array_length(ids,1) + LOOP + return QUERY + select cdb_crankshaft.CDB_DotDensity(geoms[i], coalesce(perGeom[i],0)::INTEGER); + END LOOP; +END $$ -LANGUAGE plpgsql VOLATILE; +LANGUAGE plpgsql; + + +-- +-- Daysymetric Dot Density +-- +-- @param geom: the geometry that has the +-- +-- @param no_points: the total number of points to create +-- +-- @param targetGeoms: the geometry that has the +-- +-- @param weights: targetGeom weights +-- +-- RETURNS setof points + +CREATE OR REPLACE FUNCTION CDB_DasymetricDotDensity(geom GEOMETRY, no_points NUMERIC, targetGeoms GEOMETRY[], weights numeric []) +RETURNS setof GEOMETRY +AS $$ + BEGIN + RETURN QUERY + SELECT cdb_crankshaft._CDB_WeightedDD(no_points, array_agg( ST_INTERSECTION(geom,g)), array_agg(ST_AREA(ST_INTERSECTION(geom,g))*w)::NUMERIC[]) + FROM unnest(targetGeoms) as g , unnest(weights) as w + WHERE geom && g; + END +$$ +LANGUAGE plpgsql;