-- If you are foolish enough to want to redo the above project at the
-- tract level (for california alone) this will get you
-- started. Herein we: (1) create a spatially aware table of census
-- tract polygons (2) affix to it a variable caled "logrecno" which
-- coresponds to a like named variable in each of the 76 huge files of
-- census tract level data in ICPSR_13346. ICPSR_13346 are the
-- slightly repackaged tract level data for california from the 200
-- Census.

-- This of course is only the first step. With a spatially aware table
-- of census tract polygons *with* a convenient logrecno variable in
-- place, you *still* must sift through the 76 files in ICPSR_13346;
-- read one or more into postgresql and then create the cartogram grid
-- therefrom.  This file includes a trick for extracting a selected
-- few columns froma .CSV file -- a useful trick for this and other
-- projects.

-- STEP 1 (The result of this step is in demog.215.public so you don't
-- need to do this yourself)

--- creating a table of california census tracts 


-- shp2pgsql -s 4269 -c /data/commons/carlm/USA_Counties_shapefile/tl_2009_us_county public.usa_tiger_counties >
--   /tmp/usa_county.sql

shp2pgsql -s 4269 -c /data/nobackup/carlm/06_CALIFORNIA/tl_2009_06_tract00 public.ca_tiger_tracts >  /tmp/ca_tracts.sql 

-- From the psql prompt :
set client_encoding to 'LATIN1';
\i /tmp/usa_county.sql

-- The above commands set the client encoding to match what the US gov
-- did in these shapefiles. (odd move on their part) and then reads in
-- the sql file.


-- Once read in, we now change add a new geom column where the data
-- are in a more useful and familiar SRID

-- add a new geo column making sure it gets registered in postgis
-- NOTICE that this is going into the public schema YOU should not
-- have permission to add things to the public schema so if you are
-- doin this for yourself change the scheme ok

select AddGeometryColumn('public', 'ca_tiger_tracts','geom',2163,'MULTIPOLYGON',2);

-- populate the new column with reprojection of the old the_geom column
update public.ca_tiger_tracts
 set geom= ST_Transform(the_geom,2163);

-- verify the existance in the table of both columns
\d ca_tiger_tracts

-- drop the old column
select DropGeometryColumn('public','ca_tiger_tracts', 'the_geom');



-- Now another stupid little complication: counties in CA and
-- elsewhere have islands generally these are small insignificnt
-- places. BUT because our data is by county they will cause large
-- significant problems if they are left the geom.

-- first a table with one row for each polygon in usa_tiger_counties
-- in other words we separate the MULTIPOLYGONS that are counties into
-- multiple rows containing a single plain old POLYGON

select gid,
ST_GeometryN(
	geom,
	generate_series(1,ST_NumGeometries(geom))) as geom
into temp allpoly
from public.ca_tiger_tracts;

-- how many extra bits/ islands?
select (select count(*) from allpoly) -
       (select count(*) from public.ca_tiger_tracts) ;

-- add a column holding the area in order to simplify the SQL
alter table allpoly add column area float;

update allpoly
set area = ST_Area(geom);

-- create a table that holds only the gid and the geom of the larges
-- polygon in the set of polygons that defines each county
select allpoly.* 
into temp maxpoly
from allpoly
     where area = (select max(area)
                            from allpoly as ap where ap.gid = allpoly.gid);

-- create a new version of ca_tiger_tracts  that has all of the latter
-- plus a geom1 column which holds the polygon that defines the largest bit

select ca.* , mx.geom as geom1 
into temp ca_tiger_tracts2
from public.ca_tiger_tracts as ca
inner join maxpoly as mx on ca.gid = mx.gid;

-- Let's see which counties have lost their smaller bits
select name from usa_tiger_counties2
where  not geom = geom1;


-- replace usa_tiger_counties with usa_tiger_counties2
update  ca_tiger_tracts2
set geom = geom1;

alter table ca_tiger_tracts drop column geom1;

drop table public.ca_tiger_tracts;

select * into public. ca_tiger_tracts
from ca_tiger_tracts2;

drop table  ca_tiger_tracts2;


create INDEX idx_ca_tiger_tracts ON public.ca_tiger_tracts USING gist(geom);

SELECT Populate_Geometry_Columns();

-- STEP 2 Create a table that contains enough information to match
-- with ca_tiger_tracts **AND** also contains that "logrecno" variable
-- that we'll need in order to acces the tract level data that is
-- hidden in 76 huge files


create table public.ca_tract_geo (
 logrecno varchar(20),
 sumlevel varchar(8),
 geocomp varchar(11),
 state varchar(2),	
 county varchar(3),
 tract varchar(6))

-- run this perl program on DS0077/13346-0077-Data-Census_2000_Geo._Header.txt
-- redirecting output into /tmp/GeoCols.txt

#! /usr/bin/perl -w
##Mon Nov 29 16:53:09 PST 2010
## This will read the geo summary file and produce a more manageable
## file coma separated values and containing State,county,tract and logrecid
########################################################################

use strict;

foreach(<STDIN>){
  my $state=substr($_, 6, 2);
  my $sumlev=substr($_,8, 3);
  my $geocomp=substr($_,11, 2);
  my $chariter=substr($_,13,3);
  my $county=substr($_,31,3);
  my $tract=substr($_,55,6);
  my $logrecno=substr($_,18,7);

  if($tract !~ /^\s+$/){
    ## tract is not blank
    print(join(',',$logrecno,$sumlev,$geocomp,$state,$county,$tract),"\n");
  }
}
--- End of perl program

\copy public.ca_tract_geo from /tmp/GeoCols.txt DELIMITER ','

create index idx_ca_tract_geo1 on public.ca_tract_geo(logrecno);
create index idx_ca_tract_geo2 on public.ca_tract_geo(state);
create index idx_ca_tract_geo3 on public.ca_tract_geo(county);

-- add a useful column to ca_tiger_tracts
alter table public.ca_tiger_tracts
 add column logrecno varchar(20);

-- sumlevel is the "sumaraly level" varaible which the census uses
-- with each record to describe the geography.  It is quite confusing
-- because the geopgraphies overlap and so some codes indicate partial
-- coverage. Don't worry too much about it.  the value "140" is the
-- one that says complete census tract.

update public.ca_tiger_tracts tig
 set logrecno = (select logrecno from public.ca_tract_geo geo where
                 geo.tract = tig.tractce00 AND
		 geo.county = tig.countyfp00 AND
		 geo.sumlevel = '140');

 create index idx_ca_tiger1 on public.ca_tiger_tracts(logrecno);

-- So now what?

-- California tract level data is in ~carlm/215/GIS/ICPSR_13346. In
-- that directory you will find a data dictionary and 77
-- subdirectories.  Each subdirectory holds one segment of data DS001
-- for example holds. Luckily for you, other than the geo header data
-- in DS0077, which I have dealt with, the data are all in CSV
-- files. This makes it easy to extract just the columns that you are
-- interested in OTHERWISE you would need to write a create statement
-- that included gobs of variables that you do not care about.

-- HERE IS THE TRICK: suppose that you only care about colum 5, 9,and
-- 22 of the file in DS0002.  The file
-- ~carlm/215/GIS/ICPSR_13346/13346-Data_ditionary.txt holds the
-- secret of what those columns mean: 5 = logrecno; 9= number of
-- married-couple families with children under 6 only.  22= number of
-- female householder-no husband families with children under 6 only.

-- Extract the 3 interesting columns into a /tmp file
-- In a terminal window:
cd ~carlm/215/GIS/ICPSR_13346/DS0002

awk -F','  '{print $5,$9,$22}' OFS=',' \
   <13346-0002-Data-Census_2000_File_Seg._2.txt > \
   /tmp/selectedTractData.csv

-- Read the file from /tmp into postgresql
create table tract_data0 (
  logrecno varchar(20),
  mcouples int,
  femheads int);

\copy  tract_data0 from /tmp/selectedTractData.csv DELIMITER ','

select count(*) from ca_tiger_tracts;

select count(*) from ca_tiger_tracts tig, tract_data0 t0
where tig.logrecno = t0.logrecno;

-- Once you have a table with the variable(s) that you would like to
-- use for "densisty" in your cartogram, you can pick up the thread at
-- http://www.demog.berkeley.edu/215/GIS/02.1CartIN.sql which is the
-- solid box in Step 2: Creating a grid of points ...