
-- This file contains all the tedious steps involved in creating a postGIS
-- table that has ONE polygon for each county (no islands). The table 
-- is in the pubic schema of demog215

-- The necessary steps are:
-- 1) find and dowload the shapefile from the Census Bureau

-- 2) read the shape file into PostGIS

-- 3) convert the geographic data into the SRID that we have been using

-- 4) convert the MULTIPOLYGONS into simple polygons -- some counties
-- have islands and those will screw up that cartogram royally unless
-- they are removed.

-- 5) create indexes and populate the geometry_columns table so that
-- we can run queries fast and possibly use the table in OpenJump
-- et. al.


-- STEP 1
-- I downloaded this zip file
-- ftp://ftp2.census.gov/geo/tiger/TIGER2009/tl_2009_us_county.zip
-- Tiger files are stored in SRID=4269 which is almost the same lat
-- and long -- in It will be easier for us to work with the counties
-- if they are stored in the same projection/SRID as the fastfoods
-- data so we will need to read the data into postGIS and then convert
-- it to the proper SRID


-- STEP 2
-- The data arrive as "shape" files.  This is the most likely format
-- that your will find polygon data in so it is convenient that
-- postGIS comes with a tool for reading such data: shp2psql 

-- at the Unix prompt:

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

-- The above line creates a file of SQL commands, in /tmp/usa_county.sql  that will create a
-- table in postgres that holds the data stored in the shaefile(s).

-- 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.

-- STEP 3
-- 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', 'usa_tiger_counties','geom',2163,'MULTIPOLYGON',2);

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

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

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


-- STEP 4
-- 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.usa_tiger_counties;

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

-- 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 usa_tiger_counties  that has all of the latter
-- plus a geom1 column which holds the polygon that defines the largest bit
select usa.* , mx.geom as geom1 
into temp usa_tiger_counties2
from public.usa_tiger_counties as usa
inner join maxpoly as mx on usa.gid = mx.gid;

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

-- Los Angeles : geom1 should be smaller than geom
select ST_Area(geom), ST_Area(geom1)
from usa_tiger_counties2 where name = 'Los Angeles';

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

alter table usa_tiger_counties2 drop column geom1;

drop table public.usa_tiger_counties;

select * into public. usa_tiger_counties
from usa_tiger_counties2;

drop table  usa_tiger_counties2;

-- STEP 5
create INDEX idx_usa_tiger_counties ON public.usa_tiger_counties USING gist(geom);

SELECT Populate_Geometry_Columns();

