# This python script is designed to combine data from the Census
# block files with data from the USGS.  The combined shape files 
# will then be read in by stata.   

# Import arcpy module
import arcpy
import shutil
import os
import arcgisscripting
gp = arcgisscripting.create()

# define datasets.  You will need to change the names and 
# directories here.  

block_alameda = "V:/atalay/gis/tl_2010_06001_tabblock10.shp"
block_contra_costa = "V:/atalay/gis/tl_2010_06013_tabblock10.shp"
block_marin = "V:/atalay/gis/tl_2010_06041_tabblock10.shp"
block_san_francisco = "V:/atalay/gis/tl_2010_06075_tabblock10.shp"
block_san_mateo = "V:/atalay/gis/tl_2010_06081_tabblock10.shp"
elev_contour = "V:/atalay/gis/elev_contour.shp"

# names of the shape files that we want to create

intersectAlameda="V:/atalay/gis/intersect_alameda.shp"
intersectContraCosta="V:/atalay/gis/intersect_contra_costa.shp"
intersectMarin="V:/atalay/gis/intersect_marin.shp"
intersectSanFrancisco="V:/atalay/gis/intersect_san_francisco.shp"
intersectSanMateo="V:/atalay/gis/intersect_san_mateo.shp"

# check to see if the file already exists, if so delete it

if gp.Exists(intersectAlameda): arcpy.Delete_management(intersectAlameda)
if gp.Exists(intersectContraCosta): arcpy.Delete_management(intersectContraCosta)
if gp.Exists(intersectMarin): arcpy.Delete_management(intersectMarin)
if gp.Exists(intersectSanFrancisco): arcpy.Delete_management(intersectSanFrancisco)
if gp.Exists(intersectSanMateo): arcpy.Delete_management(intersectSanMateo)

# produce the shape files, combining data on the 
# block groups from the Census website with the USGS data.

arcpy.Intersect_analysis([block_alameda, elev_contour], intersectAlameda, "ALL", ".01 Kilometers")
arcpy.Intersect_analysis([block_contra_costa, elev_contour], intersectContraCosta, "ALL", ".01 Kilometers")
arcpy.Intersect_analysis([block_marin, elev_contour], intersectMarin, "ALL", ".01 Kilometers")
arcpy.Intersect_analysis([block_san_francisco, elev_contour], intersectSanFrancisco, "ALL", ".01 Kilometers")
arcpy.Intersect_analysis([block_san_mateo, elev_contour], intersectSanMateo, "ALL", ".01 Kilometers")