# Script to isolate specific analysis cities

# Michael Minn - 31 May 2014

library("rgdal")
library("maptools")

cities = list(
	list(layer = "tl_2013_04_place", city = "Phoenix"),
	list(layer = "tl_2013_06_place", city = "Los Angeles"),
	list(layer = "tl_2013_06_place", city = "San Diego"),
	list(layer = "tl_2013_17_place", city = "Chicago"),
	list(layer = "tl_2013_36_place", city = "New York"),
	list(layer = "tl_2013_42_place", city = "Philadelphia"),
	list(layer = "tl_2013_48_place", city = "Houston"),
	list(layer = "tl_2013_48_place", city = "San Antonio"))

outplaces = NULL

for (city in cities)
{
	print(city$city)
	places = readOGR(dsn=".", layer=city$layer)
	place = places[places$NAME == city$city,]
	if (is.null(outplaces))
		outplaces = place
	else
		outplaces = spRbind(outplaces, place)
}

writeOGR(outplaces, dsn=".", layer="cities", driver="ESRI Shapefile", overwrite_layer = T)
