
-- Now on to creating the input grid for the cartogram:

-- First we must create a "grid" of evenly spaced points across the
-- US. Then each needs to have some value such as population assigned to
-- it.  That grid then get's imported into R to have the diffusion
-- cartogram program run upon it.

-- STEP 1
-- creating the grid of evenly spaced points

-- What we actually create below is a table with one row for each
-- combination of X and Y (lattitude and longitude if you like...but
-- actually X and Y coordinates of the SRID system that our geom is
-- stored in) that will comprise our grid.  We can decide how fine a
-- grid to use for now let's have make the grid 1025X513  since the US
-- (minus HI, AK, guam, Saipan, Puerto Rico & etc) is wider than it is
-- tall

-- first min and max extents of x and y and step size

select  (max(st_xmax(geom))::int) as Xmx,
        (min(st_xmin(geom))::int)as Xmn, 
        (max(st_ymax(geom))::int) as Ymx,
        (min(st_ymin(geom))::int) as Ymn,
 (max(st_xmax(geom))::int - min(st_xmin(geom))::int)/1025 as Xstep,
 (max(st_ymax(geom))::int - min(st_ymin(geom))::int)/513 as Ystep
into stepsize
from usa_tiger_counties where statefp::int <= 56 and statefp::int != 2
and statefp::int != 15;
 
--  statefp 15 is HI, 2 is AL above 56 are places like Guam and Puerto Rico


-- create the grid like table
select xtab.* , ytab.* into us_grid
from 
(select generate_series(
        (select Xmn from stepsize),
        (select Xmx from stepsize),
	(select Xstep from stepsize)) as xcord) as xtab,

(select generate_series(
        (select Ymn from stepsize),
        (select Ymx from stepsize),
	(select Ystep from stepsize))as ycord) as ytab ;


-- add a geom so we can determine which county each point is in
select AddGeometryColumn(SCHEMA, 'us_grid','geom',2163,'POINT',2);

-- Note that xcord refers to longitude, one crosses longitude lines as one moves in the x direction
update us_grid
set geom= ST_GeomFromText('POINT(' || xcord || ' ' || ycord || ')', 2163);

CREATE INDEX idx_us_grid_geom ON us_grid USING gist(geom);

vacuum analyze ;

-- Now let's add at least one column with which we might perform the
-- cartogram dance

-- In the public schema is a table called dc_2009_est This table
-- contains 2009 census estimates of about 200 census values for each
-- county.  We now join that data to usa_tiger_counties to create a
-- large TEMPORARY table from which we will soon extract the "density"
-- measture for each cell of the evenly spaced grid that we just
-- contructed.

----------------------------------------------------------------------
----------------------------------------------------------------------
-- NOTE:   WHILE WE CAN ALWAYS USE total population density... There
-- might well exist something clever.  For example if you have data on
-- sales of acne medication you might want to run the cartogram
-- program over population density of adolescents.  What you decide to
-- use should be a densisty measure -- so make sure it gets divided by
-- the area of the county.

-- For the present fastfood restaurant question perhaps you might wish
-- to add some alternative density measures...?
----------------------------------------------------------------------
----------------------------------------------------------------------

\d dc_2009_est

-- hmmm what might those columns mean?

select * from dc_vnames;


-- For want of a clever name.. we now create "bigfile" a file that
-- will hold both the tiger geography data and census estimates for 2009.

select  tig.*, dc.* ,ST_Area(tig.geom) as area 
into temp bigfile
from public.usa_tiger_counties as tig, public.dc_2009_est as dc 
where tig.cntyidfp = dc.geo_id2;

-- how many records are missing from bigfile?
select  (select count(*) from bigfile )-
        (select count(*) from usa_tiger_counties);

-- and which ones are they?
SELECT cntyidfp, name, name from usa_tiger_counties
where cntyidfp  not in (select geo_id2 from dc_2009_est);
 
-- guam,puerto rico, saipan, ...


-- Now back to the us_grid: we would like to put a number in each row
-- of us_grid that coresponds to some "density" measure (total
-- population for example) of the county in which the point lies. This
-- is the key input to the cartogram procedure.

CREATE INDEX idx_bigfile_geom ON bigfile USING gist(geom);

vacuum analyze ;

-- STEP 2
-- create the cart_in table each row of which holds an element in the
-- evenly spaced grid that we talked of earlier AND some values that
-- we might want to use  as "density" measures.

-- even though it generates a warning we use a left innner join b/c we
-- want a complete grid even though many points in us_grid will fall
-- into the Atlantic, Pacific, Mexico or Canada.

-- ST_Intersects tests for any common point in two geometries. What we
-- want to test for in this query is is whether a point in us_grid is
-- inside of a polygon of a particular county.  There are other ST_
-- functions that will do this

 --ST_Covers( bf.geom,ug.geom); or  ST_Within(ug.geom,bf.geom);

-- but ST_Intersects turns out to be the fastest when one is looking
-- to verify that a point is within a polygon.
 
-- This one takes more than a minute to run.... 
select ug.*, bf.statefp, bf.countyfp, bf.name,bf.gid, bf.area,
 bf.t001001 as totpop, (bf.t001001/bf.area)::float as totpopden,
 bf.T0082009003 as boys , (bf.T0082009003/bf.area)::float as boysdens 
into  cart_in
    from  us_grid as ug INNER JOIN 
    	     bigfile as bf on
	     ST_Intersects (ug.geom, bf.geom);

CREATE INDEX idx_cart_in_geom ON cart_in USING gist(geom);


SELECT Populate_Geometry_Columns();

