Create a random point inside a Census block, one per person

random_points_census_blocks.sql

-- Function: public.pop_tract_explode()
-- DROP FUNCTION public.pop_tract_explode();

CREATE OR REPLACE FUNCTION public.pop_block_explode()
  RETURNS text AS
$BODY$
DECLARE
    record    RECORD;

    gid INTEGER;
    geom_2232 GEOMETRY;
    geoid10 BIGINT;

    -- Variables for explode
    population INTEGER;
	pop_white INTEGER;
	pop_black INTEGER;
	population_divided INTEGER;
	divider INTEGER := 1;

    category_tag TEXT;   
    ct integer;
    skipped INTEGER := 0;

    -- Variables for randomize steps
    i INTEGER;
    maxiter INTEGER := 1000000;
    x0 DOUBLE PRECISION;
    dx DOUBLE PRECISION; 
    y0 DOUBLE PRECISION;
    dy DOUBLE PRECISION;
    xp DOUBLE PRECISION;
    yp DOUBLE PRECISION;
    poly_geom GEOMETRY;
    randPoint GEOMETRY;
	--randPoint_3435 GEOMETRY;

BEGIN       

    -- Creating a destination table with a matching data model to the input table and adding needed fields
    EXECUTE 'DROP TABLE IF EXISTS pop_block_explode';
    EXECUTE 'CREATE TABLE pop_block_explode AS 
             SELECT "gid", "geoid10", "pop_total", "pop_white", "pop_black" FROM public.chicago_block_with_population LIMIT 1;';
    EXECUTE 'TRUNCATE TABLE pop_block_explode;';
	EXECUTE 'ALTER TABLE pop_block_explode ADD COLUMN "population_divided" NUMERIC;';
	EXECUTE 'ALTER TABLE pop_block_explode ADD COLUMN "divider" NUMERIC;';
    --EXECUTE 'ALTER TABLE pop_block_explode ADD COLUMN "category_tag" character varying(100);';
    EXECUTE 'ALTER TABLE pop_block_explode ADD COLUMN rand_point GEOMETRY;';
	--EXECUTE 'ALTER TABLE pop_block_explode ADD COLUMN rand_point_3435 GEOMETRY;';

    -- Loop through all records of the input table
    FOR record IN EXECUTE 'SELECT gid, geoid10, ST_Transform(geom, 3435) AS geom, "pop_total", "pop_white", "pop_black" FROM public.chicago_block_with_population;'
	LOOP	
		-- Set the polygon geometry variable for the current record
		poly_geom := record.geom;

		-- Determine the bounds of polygon geometry
		x0 = ST_XMin(poly_geom);
		dx = (ST_XMax(poly_geom) - x0);
		y0 = ST_YMin(poly_geom);
		dy = (ST_YMax(poly_geom) - y0);	

		-- Set the variables for the current record
		gid := record.gid::INTEGER;
		geoid10 := record.geoid10::BIGINT;	
		population := record.pop_total::INTEGER;
		pop_white := record.pop_white::INTEGER;
		pop_black := record.pop_black::INTEGER;
		--population_divided := round(record.population::INTEGER/divider);

		-- Loop throu
		ct := 0;
		category_tag := 'population';
		WHILE ct < population::decimal/divider  -- you can make each point represent a specific number of occurances by dividing by that value (i.e. population/100 makes 1 pt = 100 people)
		LOOP
			RAISE NOTICE 'Count: %', ct;

			-- Only try to place the point within the polygon if not exceeding the maximum interation value
			i := 0;		
			WHILE i < maxiter 
			LOOP
				i = i + 1;
				xp = x0 + dx * random();
				yp = y0 + dy * random();
				randPoint = ST_SetSRID( ST_MakePoint( xp, yp ), ST_SRID(poly_geom) );
				--randPoint_3435 = ST_SetSRID( ST_Transform(ST_MakePoint( xp, yp ), 3435), ST_SRID(poly_geom) );

				EXIT WHEN ST_Within( randPoint, poly_geom ); 
			END LOOP;

		 	IF i >= maxiter THEN
				RAISE NOTICE 'RandomPoint: number of interations exceeded %', maxiter;
				skipped = skipped + 1;
			ELSE
			-- Insert a new record into the destination table for the random point 					
			INSERT INTO pop_block_explode VALUES ( gid,geoid10,population,pop_white,pop_black,round(population::decimal/divider),divider,randPoint );
			END IF;

			ct = ct + 1;
		END LOOP;
	END LOOP;

	-- Create spatial index on the new spatial column
	EXECUTE 'CREATE INDEX pop_block_explode_gist ON pop_block_explode USING gist (rand_point)';

    RETURN skipped;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION public.pop_block_explode()
  OWNER TO stevevance;