classify a list of geometries via table; extract grid cells for each class; build indexes and add pkeys

do_make_region_grids.py

tSQL = '''
drop table if exists region_{0}_grid cascade;
create table region_{0}_grid as
select distinct on (g.gridcode) g.gridcode as gid, 
  st_transform( ST_GeometryN(g.the_geom,1), 4326)::geometry(Polygon,4326) as geom 
from grid150m_uniq g, 
  tl_2010_06_county10_3310 c, region_defs r
where st_intersects( g.the_geom, c.geom ) AND 
  c.countyfp10 = r.fips AND
  r.region_id = {0};

alter table region_{0}_grid add primary key(gid);
create index r{0}g_geom_idx on region_{0}_grid using GIST (geom);
'''

tSQL_3310 = '''
drop table if exists region_{0}_grid_3310 cascade;
create table region_{0}_grid_3310 as
select
  o.gid,
  ST_transform(o.geom, 3310)::geometry(Polygon,3310) as geom
from region_{0}_grid o;

alter table region_{0}_grid_3310 add primary key(gid);
create index r{0}g3_geom_idx on region_{0}_grid_3310 using GIST (geom);
'''

import psycopg2

try:
  conn = psycopg2.connect( "dbname=grid_load" )
  curs = conn.cursor()
except Exception, E:
  print str(E)
  exit(1)

for n in xrange(1,9):
  try:
    curs.execute( tSQL.format(n) )
  except Exception, E:
    print str(E)
  conn.commit()

for n in xrange(1,9):
  try:
    curs.execute( tSQL_3310.format(n) )
  except Exception, E:
    print str(E)
  conn.commit()