library("rgdal")
# library("rgeos")

wgs84 = CRS("+init=epsg:4326")

create_census_zip_table = function()
{
	# ZIP Codes
	if (!exists("zip_codes"))
		zip_codes <<- readOGR(dsn="/home/michael/mapping-sustainability/census-zip/downloads", 
			layer="tl_2013_us_zcta510")
		zip_codes = spTransform(zip_codes, CRS=wgs84)

	centroids <<- SpatialPoints(coordinates(zip_codes), proj4string=wgs84)


	# Find County Names

	if (!exists("counties"))
		counties <<- readOGR(dsn="/home/michael/mapping-sustainability/boundaries/analysis-areas",
			layer="boundaries-county")
		counties = spTransform(counties, CRS=wgs84)

	analysis_polygons <<- SpatialPolygons(counties@polygons, proj4string=wgs84)

	indices <<- over(centroids, analysis_polygons)

	zip_codes@data$COUNTY <<- as.character(counties@data[indices, "NAMELSAD10"])


	# Find City Names

		cities <<- readOGR(dsn="/home/michael/mapping-sustainability/boundaries/analysis-areas",
			layer="boundaries-city")
		cities = spTransform(cities, CRS=wgs84)

	analysis_polygons <<- SpatialPolygons(cities@polygons, proj4string=wgs84)

	indices <<- over(centroids, analysis_polygons)

	zip_codes@data$CITY <<- as.character(cities@data[indices, "NAME"])


	# Keep only zip codes within city or county analysis areas
	census_zip <<- zip_codes[((!is.na(zip_codes@data$CITY)) | (!is.na(zip_codes@data$COUNTY))), ]

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

add_acs_columns = function()
{
	# If just the data slots are merged, the polygon/attribute order gets messed up
	# x@data = merge(x@data, y) = bad

	# DP03 - SELECTED ECONOMIC CHARACTERISTICS: 2008-2012 American Community Survey 5-Year Estimates
	dp03_columns = c(
		"HC03_VC13" = "PCUNEMPLOY",
		"HC03_VC29" = "PCDRIVE",
		"HC01_VC85" = "HHINCOME")

	print("Merging DP03")
	dp03 = read.csv("downloads/ACS_12_5YR_DP03.csv", fill=T)
	dp03 = dp03[,c("GEO.id2", names(dp03_columns))]
	colnames(dp03) = c("GEO.id2", dp03_columns)
	census_zip <<- merge(x=census_zip, dp03, by.x = "GEOID10", by.y = "GEO.id2")

	# DP04 - SELECTED HOUSING CHARACTERISTICS: 2008-2012 American Community Survey 5-Year Estimates
	dp04_columns = c(
		"HC03_VC64" = "HHRENTER",
		"HC03_VC82" = "PCNOVEHICL")

	print("Merging DP04")
	dp04 = read.csv("downloads/ACS_12_5YR_DP04_with_ann.csv", fill=T)
	dp04 = dp04[-1,] # delete annotation row
	dp04 = dp04[,c("GEO.id2", names(dp04_columns))]
	for (col in names(dp04_columns)) dp04[,col] = as.numeric(as.character(dp04[,col]))
	colnames(dp04) = c("GEO.id2", dp04_columns)
	census_zip <<- merge(x=census_zip, dp04, by.x = "GEOID10", by.y = "GEO.id2")


	# DP05 - ACS DEMOGRAPHIC AND HOUSING ESTIMATES: 2008-2012 American Community Survey 5-Year Estimates 
	# "Race alone or in combination with one or more other races"
	dp05_columns = c(
		"HC01_VC03" = "POPTOTAL",
		"HC01_VC21" = "MEDIANAGE",
		"HC03_VC72" = "PCWHITE", 
		"HC03_VC73" = "PCBLACK",
		"HC03_VC82" = "PCHISPANIC")

	print("Merging DP05")
	dp05 = read.csv("downloads/ACS_12_5YR_DP05_with_ann.csv", fill=T)
	dp05 = dp05[-1,] # delete annotation row
	dp05 = dp05[,c("GEO.id2", names(dp05_columns))]
	for (col in names(dp05_columns)) dp05[,col] = as.numeric(as.character(dp05[,col]))
	colnames(dp05) = c("GEO.id2", dp05_columns)
	census_zip <<- merge(x=census_zip, dp05, by.x = "GEOID10", by.y = "GEO.id2")

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

add_leed_location_count = function()
{
	print("Merging LEED Project location counts")
	locations <<- read.csv("/home/michael/mapping-sustainability/green-buildings/leed-locations.csv", fill=T)
	locations$postal_code = substr(locations$postal_code, 1, 5) # clip any zip+4

	codes = as.data.frame(table(locations$postal_code))
	colnames(codes) = c("ZIP", "LEED")
	census_zip <<- merge(x = census_zip, y = codes, by.x = "GEOID10", by.y = "ZIP")

	codes = as.data.frame(table(locations[locations$brownfield == "yes","postal_code"]))
	colnames(codes) = c("ZIP", "LEEDBROWN")
	census_zip <<- merge(x = census_zip, y = codes, by.x = "GEOID10", by.y = "ZIP")

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

add_bike_rack_count = function()
{
	print("Adding Bike Racks")

	# Load data
	racks <<- readOGR(dsn="/home/michael/mapping-sustainability/bike-sharing", layer="bike-sharing")

	# racks@data$BIKERACKS <<- NULL

	# Find ZIP code for each rack
	rack_points <<- SpatialPoints(coordinates(racks), proj4string=CRS(proj4string(racks)))
	zip_polygons <<- SpatialPolygons(census_zip@polygons, proj4string=CRS(proj4string(census_zip)))
	zip_polygons <<- spTransform(zip_polygons, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
	indices <<- over(rack_points, zip_polygons)
	racks@data$ZIP <<- as.character(census_zip@data[indices, "GEOID10"])


	# Aggregate to find counts by zip code
	sites = aggregate(Docks ~ ZIP, data = racks@data, FUN = length)
	docks = aggregate(Docks ~ ZIP, data = racks@data, FUN = sum)

	# Combine aggregations and give meaningful shapefile column names
	colnames(sites) = c("ZIP", "BIKESITES")
	colnames(docks) = c("ZIP", "BIKERACKS")

	# Merge with zip code shapefile
	census_zip <<- merge(x = census_zip, y = sites, by.x = "GEOID10", by.y = "ZIP")
	census_zip <<- merge(x = census_zip, y = docks, by.x = "GEOID10", by.y = "ZIP")

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

add_rain_barrel_count = function()
{
	print("Adding rain barrels")

	# Load data
	csv = read.csv("/home/michael/mapping-sustainability/rain-barrels/chicago-rain-barrels.csv", fill=T, check.names=F)
	csv = csv[,c("Rebate Applicant: Zip Code","Barrel Capacity (gallons)")]
	colnames(csv) = c("ZIP", "GALLONS")

	# Aggregate to find total count and gallons by zip code
	count = aggregate(GALLONS ~ ZIP, data = csv, FUN = length)
	gallons = aggregate(GALLONS ~ ZIP, data = csv, FUN = sum)

	# Load csv
	csv = read.csv("/home/michael/mapping-sustainability/rain-barrels/san-diego-rain-barrels.csv", fill=T, check.names=F)
	csv = csv[,c("Zip Code","Capacity (gallons)")]
	colnames(csv) = c("ZIP", "GALLONS")

	# Aggregate to find total count and gallons by zip code
	count2 = aggregate(GALLONS ~ ZIP, data = csv, FUN = length)
	gallons2 = aggregate(GALLONS ~ ZIP, data = csv, FUN = sum)

	# Load csv and overlay to find ZIP code
	csv = read.csv("/home/michael/mapping-sustainability/rain-barrels/philadelphia-RainBarrels_130923.csv",
		fill=T, check.names=F)
	csv_points = SpatialPoints(csv, proj4string=CRS(proj4string(census_zip)))
	zip_polygons = SpatialPolygons(census_zip@polygons, proj4string=CRS(proj4string(census_zip)))
	indices = over(csv_points, zip_polygons)
	csv$ZIP = as.character(census_zip@data[indices, "GEOID10"])

	# Aggregate to find total count and gallons by zip code
	count3 = as.data.frame(table(csv$ZIP))
	colnames(count3) = c("ZIP", "GALLONS")
	gallons3 = count3
	gallons3$GALLONS = NA

	# Combine aggregations and give meaningful shapefile column names
	count = rbind(count, count2, count3)
	gallons = rbind(gallons, gallons2, count3)
	colnames(count) = c("ZIP", "RAINBARREL")
	colnames(gallons) = c("ZIP", "RAINGALLON")

	# Merge with zip code shapefile
	# census_zip@data$RAINBARREL = NULL
	census_zip <<- merge(x = census_zip, y = count, by.x = "GEOID10", by.y = "ZIP")

	# Gallons are proportional to count since barrels are generally the same size
	# census_zip <<- merge(x = census_zip, y = gallons, by.x = "GEOID10", by.y = "ZIP")
	
	#census_zip@data[is.na(census_zip@data$RAINBARREL), "RAINBARREL"] = 0
	#spplot(obj=census_zip[census_zip@data$CITY == "Chicago",], zcol="RAINBARREL")

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

add_brownfield_counts = function()
{
	layer_names = c(
		"RCRA" = "rcra",
		"CERCLIS" = "cerclis",
		"CIMC" = "cimc",
		"BROWNEPA" = "epa-grants",
		"BROWNSTATE" = "state-programs",
		"BROWNCITY" = "city-programs")

	zip_polygons <<- SpatialPolygons(census_zip@polygons, proj4string=CRS(proj4string(census_zip)))
	zip_polygons <<- spTransform(zip_polygons, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

	for (column in names(layer_names))
	{
		print(column)

		brownfields = readOGR(dsn="/home/michael/mapping-sustainability/brownfields", layer=layer_names[column])
		brownpoints = SpatialPoints(coordinates(brownfields), 
			proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

		indices = over(brownpoints, zip_polygons)
		brownzips = as.character(census_zip@data[indices, "GEOID10"])

		browntable = as.data.frame(table(brownzips))
		colnames(browntable) = c("ZIP", column)
		census_zip <<- merge(x = census_zip, y = browntable, by.x = "GEOID10", by.y = "ZIP")
	}

	writeOGR(census_zip, dsn="/home/michael/mapping-sustainability/census-zip", 
		layer="census-zip", driver="ESRI Shapefile", overwrite_layer=T)
}

create_census_zip_table()

add_acs_columns()

add_leed_location_count()

# census_zip <<- readOGR(dsn="/home/michael/mapping-sustainability/census-zip", layer="census-zip")

add_bike_rack_count()

add_rain_barrel_count()

add_brownfield_counts()

