library(rgdal)
library(maptools)
library(spatstat)
library(raster)

# Load data

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

if (!exists("analysis_area"))
{
	analysis_area = readOGR(dsn="/home/michael/mapping-sustainability/boundaries/analysis-areas", layer="boundaries-city")
	analysis_area = analysis_area[(analysis_area$NAME == "Chicago"), ]
	analysis_area = spTransform(analysis_area, CRS=wgs84)
	analysis_area = SpatialPolygons(analysis_area@polygons, proj4string=wgs84)

	width = 512
	height = width * (analysis_area@bbox["y","max"] - analysis_area@bbox["y", "min"]) /
		(analysis_area@bbox["x","max"] - analysis_area@bbox["x", "min"])

	green_buildings = readOGR(dsn="/home/michael/mapping-sustainability/green-buildings", layer="leed-locations")
	green_buildings = spTransform(green_buildings, CRS=wgs84)
	indices = over(SpatialPoints(green_buildings@coords, proj4string=wgs84), analysis_area)
	green_buildings = green_buildings[!is.na(indices), ]

	green_buildings@data[is.na(green_buildings$brwnfld),"brwnfld"] = "no"
	brown_buildings = green_buildings[green_buildings$brwnfld == "yes",]

	bike_racks = readOGR(dsn="/home/michael/mapping-sustainability/bike-sharing", layer="bike-sharing")
	bike_racks = spTransform(bike_racks, CRS=wgs84)
	indices = over(SpatialPoints(bike_racks@coords, proj4string=wgs84), analysis_area)
	bike_racks = bike_racks[!is.na(indices), ]

	brownfields = readOGR(dsn="/home/michael/mapping-sustainability/brownfields", layer="state-programs")
	brownfields = spTransform(brownfields, CRS=wgs84)
	indices = over(SpatialPoints(brownfields@coords, proj4string=wgs84), analysis_area)
	brownfields = brownfields[!is.na(indices), ]

	zips = readOGR(dsn="/home/michael/mapping-sustainability/census-zip", layer="census-zip")
	zips = spTransform(zips, CRS=wgs84)
	zips@data$CITY = as.character(zips@data$CITY)
	zips[is.na(zips$CITY),"CITY"] = ""
	zips = zips[zips$CITY == "Chicago",]

	rain_barrel_zip_counts = zips$RAINBARREL
	rain_barrel_zip_counts[is.na(rain_barrel_zip_counts)] = 0
	rain_barrels = dotsInPolys(zips, f="regular", as.integer(rain_barrel_zip_counts))
	rain_barrels@proj4string = wgs84

	tracts = readOGR(dsn="/home/michael/mapping-sustainability/census-tracts", layer="tracts")
	tracts = spTransform(tracts, CRS=wgs84)
	tracts@data$CITY = as.character(tracts@data$CITY)
	tracts[is.na(tracts$CITY),"CITY"] = ""
	tracts = tracts[tracts$CITY == "Chicago",]
	tract_polygons = SpatialPolygons(tracts@polygons, proj4string=wgs84)

	hotels = read.csv(file = "hotels/chicago-yelp.csv", header=T, check.names=F)
	hotels = SpatialPointsDataFrame(coords = hotels[,c("Long","Lat")], data=hotels, proj4string=wgs84)
	indices <<- over(SpatialPoints(hotels@coords, proj4string=wgs84), analysis_area)
	hotels = hotels[!is.na(indices),]
}

# Utility function

data_table_to_html = function(table, filename = NA)
{
	if (!is.na(filename))
		sink(filename)

	cat("<html>\n<body>\n<table>\n<thead>\n<tr>\n<th>&nbsp;</th>\n")
	for (column in colnames(table))
		cat(paste(sep="", "<th>", column, "</th>\n"))

	cat("</tr>\n</thead>\n<tbody>")

	for (row in rownames(table))
	{
		cat(paste(sep="", "<tr>\n<th>", row, "</th>\n"))
		for (column in colnames(table))
			cat(paste(sep="", "<td>", table[row, column], "</td>\n"))
		cat("</tr>\n")
	}

	cat("</tbody>\n</table>\n");

	if (!is.na(filename))
		sink()
}


# Density by ZIP Code

if (0)
{
	zips$GREENBLDG = over(SpatialPolygons(zips@polygons, proj4string=wgs84), 
		SpatialPoints(green_buildings@coords, proj4string=wgs84))

	zips$GREENBROWN = over(SpatialPolygons(zips@polygons, proj4string=wgs84), 
		SpatialPoints(brown_buildings@coords, proj4string=wgs84))

	zips$HOTELS = over(SpatialPolygons(zips@polygons, proj4string=wgs84), 
		SpatialPoints(hotels@coords, proj4string=wgs84))

	census_columns = c("PCUNEMPLOY", "PCDRIVE", "HHINCOME", "HHRENTER",
		"PCNOVEHICL", "MEDIANAGE", "PCWHITE", "PCBLACK", "PCHISPANIC", "HOTELS")

	density_columns = c("GREENBLDG", "GREENBROWN", "BIKERACKS", "BROWNSTATE", "RAINBARREL")

	results = matrix(nrow = length(density_columns), ncol=length(census_columns),
			dimnames = list(density_columns, census_columns))

	for (density_column in density_columns)
		for (census_column in census_columns)
		{
			formula_string = paste(density_column, "~", census_column)
			formula = as.formula(formula_string)
			model = lm(formula, data = zips)
			print(paste(formula_string, round(summary(model)$r.squared, 3)))

			results[density_column, census_column] = sprintf("%.03f", summary(model)$r.squared)

			if (summary(model)$r.squared > 0.2)
			{
				filename = paste(sep="", "website/graphics/zip-correlation-", 
					tolower(density_column), "-", tolower(census_column), ".png")
				png(filename, width=550, height=400, res=72)

				plot(x=zips@data[,census_column], y = zips@data[,density_column], 
					xlab=census_column, ylab=density_column, frame.plot=F, pch=16, col="navy", las=1)
				abline(model, lwd=2, col="red3")
				text(x=max(zips@data[,census_column], na.rm=T) * 0.5, 
					y=max(zips@data[,density_column], na.rm=T) * 0.9,
					bquote(R^2 == .(round(summary(model)$r.squared, 3))))

				axis(side=1, lwd=3, las=1, col="white")
				axis(side=2, lwd=3, las=1, col="white")
				grid()
				dev.off()
				# readline("Enter to continue")
			}

		}

	data_table_to_html(results, "temp.html")
}


# Logit Models

if (0)
{
	tracts$HOTELS = over(tract_polygons, SpatialPoints(hotels@coords, proj4string=wgs84))

	tracts$GREENBLDG = over(tract_polygons, SpatialPoints(green_buildings@coords, proj4string=wgs84))
	tracts$GREENBLDG[is.na(tracts$GREENBLDG)] = 0
	tracts$GREENBLDG[tracts$GREENBLDG > 0] = 1

	tracts$GREENBROWN = over(tract_polygons, SpatialPoints(brown_buildings@coords, proj4string=wgs84))
	tracts$GREENBROWN[is.na(tracts$GREENBROWN)] = 0
	tracts$GREENBROWN[tracts$GREENBROWN > 0] = 1

	tracts$BIKERACKS = over(tract_polygons, SpatialPoints(bike_racks@coords, proj4string=wgs84))
	tracts$BIKERACKS[is.na(tracts$BIKERACKS)] = 0
	tracts$BIKERACKS[tracts$BIKERACKS > 0] = 1

	tracts$BROWNSTATE = over(tract_polygons, SpatialPoints(brownfields@coords, proj4string=wgs84))
	tracts$BROWNSTATE[is.na(tracts$BROWNSTATE)] = 0
	tracts$BROWNSTATE[tracts$BROWNSTATE > 0] = 1


	tracts$RAINBARREL = over(tract_polygons, SpatialPoints(rain_barrels@coords, proj4string=wgs84))
	tracts$RAINBARREL[is.na(tracts$RAINBARREL)] = 0
	tracts$RAINBARREL[tracts$RAINBARREL > 0] = 1

	census_columns = c("PCUNEMPLOY", "PCDRIVE", "HHINCOME", "HHRENTER",
		"PCNOVEHICL", "MEDIANAGE", "PCWHITE", "PCBLACK", "PCHISPANIC", "HOTELS")
	logit_columns = c("GREENBLDG", "GREENBROWN", "BIKERACKS", "BROWNSTATE", "RAINBARREL")

	results = matrix(nrow = length(logit_columns), ncol=length(census_columns), 
		dimnames = list(logit_columns, census_columns))

	for (logit_column in logit_columns)
		for (census_column in census_columns)
		{
			formula_string = paste(logit_column, "~", census_column)
			formula = as.formula(formula_string)
			model  = glm(formula, data = tracts@data, family = "binomial")
			pseudo.r.squared =  1 - (model$deviance / model$null.deviance)
			print(paste(formula_string, "Pseudo R^2 =", round(pseudo.r.squared, 3)))
			if (pseudo.r.squared >= 0.09)
			{
				print(summary(model))

				filename = paste(sep="", "website/graphics/tract-logit-", 
					tolower(logit_column), "-", tolower(census_column), ".png")
				png(filename, width=550, height=400, res=72)

				max = max(tracts@data[,census_column], na.rm=T)
				min = min(tracts@data[,census_column], na.rm=T)
				range = max - min
				breaks = min + ((0:10) * (range) / 10)

				total = table(cut(tracts@data[,census_column], breaks=breaks))
				count = table(cut(tracts@data[(tracts@data[,logit_column] > 0),census_column], breaks=breaks))
				discrete = count / total
				barplot(discrete, ylim=c(0,1), col=rgb(0, 0, 128, maxColorValue=255), border=NA,
					xlab = census_column, ylab = paste(sep="", "P(", logit_column, ")"))

				breakframe = data.frame(breaks)
				colnames(breakframe) = census_column
				logodds = predict(model, newdata=breakframe)
				odds = exp(logodds)
				probabilities = odds / (odds + 1)
				lines(probabilities, lwd=2, col= rgb(255, 204, 0, maxColorValue=255))

				text(x = 5.5, y = 0.9, labels = paste("Pseudo R^2 =", round(pseudo.r.squared, 3)))

				dev.off()
			}

			results[logit_column, census_column] = 
				paste(sprintf("%.03f", pseudo.r.squared), "/", round(model$aic, 0))
		}

	data_table_to_html(results, "temp.html")
}


# Nearest neighbor statistics between point sets

if (0)
{
	point_set_names = c("green_buildings", "brown_buildings", "bike_racks",
		"brownfields", "rain_barrels", "hotels")
	point_set_labels = c("Green Buildings", "Green Buildings on Brownfields", "Bike Racks", 
		"Remediated Brownfields", "Rain Barrels", "Hotels")

	results = matrix(ncol=length(point_set_names), nrow=length(point_set_names),
		dimnames=list(point_set_names, point_set_names))

	for (x in point_set_names)
		for (y in point_set_names)
		{
			epsg3529 = CRS("+init=epsg:5329")
			xpoints = SpatialPoints(spTransform(get(x), CRS=epsg3529))
			ypoints = SpatialPoints(spTransform(get(y), CRS=epsg3529))
			xppp = as.ppp(xpoints)
			yppp = as.ppp(ypoints)
			z = nncross(xppp, yppp)
			print(paste(x, y))
			print(summary(z))
			results[x, y] = round(mean(z$dist), 0)
			# results[x, y] = paste(round(mean(z$dist), 0), "/", round(sd(z$dist), 0))
		}

	rownames(results) = point_set_labels
	colnames(results) = point_set_labels
	data_table_to_html(results, "temp.html")
}



# Zoning Distribution

if (1)
{
	# http://secondcityzoning.org/zones
	zone_classes = c("Business", "Commercial", "Manufacturing", "Residential", "Planned Development", 
		"Planned Mfg District", "Downtown Mixed Use", "Downtown Core", 
		"Downtown Residential", "Downtown Service", "Transportation", "Parks and Open Space")

	zones = readOGR(dsn="/home/michael/mapping-sustainability/zoning", layer="Zoning_nov2012")
	zones = spTransform(zones, CRS=wgs84)
	zones$ZONE_TYPE = factor(x=zones$ZONE_TYPE, levels=1:12, labels=zone_classes)

	green_zones = over(SpatialPoints(green_buildings@coords, proj4string=wgs84), zones)

	brown_zones = over(SpatialPoints(brown_buildings@coords, proj4string=wgs84), zones)

	bike_zones = over(SpatialPoints(bike_racks@coords, proj4string=wgs84), zones)

	brownfield_zones = over(SpatialPoints(brownfields@coords, proj4string=wgs84), zones)

	barrel_zones = over(SpatialPoints(rain_barrels@coords, proj4string=wgs84), zones)

	rows = list("Green Buildings" = table(green_zones$ZONE_TYPE),
		"Green Buildings on Brownfields" = table(brown_zones$ZONE_TYPE),
		"Bike Racks" = table(bike_zones$ZONE_TYPE),
		"Remediated Brownfields" = table(brownfield_zones$ZONE_TYPE),
		"Rain Barrels" = table(barrel_zones$ZONE_TYPE))

	# png("website/graphics/zoning-distribution.png", width=550, height=400, res=72)

	par(mfrow=c(6,1))
	par(mar=c(2,4,1,2))
	for (row in names(rows))
	{
		percents = 100 * rows[[row]] / sum(rows[[row]])
		midpoints = barplot(percents, ylab=row, las=3, border=NA, xaxt = ifelse(row == "Rain Barrels", "l", "n"))
		labels = sapply(round(percents, 1), function(p) paste(p, "%", sep=""))
		text(x = midpoints, y = percents / 3, cex=0.8, pos=3, labels = labels)
	}

	# dev.off()
}


temp_tract_raster = function()
{
	tracts = readOGR(dsn="/home/michael/mapping-sustainability/census-tracts", layer="tracts")
	tracts = spTransform(tracts, CRS=wgs84)
	selection = (tracts$CITY == "Chicago")
	selection[is.na(selection)] = F
	analysis_area_tracts = tracts[selection,]

	tract_raster = raster(extent(analysis_area), ncol=width, nrow=height)
	# tract_raster = rasterize(x=analysis_area_tracts, y = tract_raster, field=analysis_area_tracts$HHINCOME / 10000)
	tract_raster = rasterize(x=analysis_area_tracts, y = tract_raster, field=analysis_area_tracts$PCDRIVE)
	tract_raster = as(tract_raster, "SpatialGridDataFrame")
	writeGDAL(tract_raster, fname="temp2.png", drivername="PNG")

	pcdrive = tract_raster@data[,1]
	bike_intensity = bike_grid@data[,1]
	plot(x=pcdrive, y=bike_intensity)
	model = lm(bike_intensity ~ pcdrive)
}

bike_rack_kernel_density = function()
{
	bike_racks = readOGR(dsn="/home/michael/mapping-sustainability/bike-sharing", layer="bike-sharing")
	bike_racks <<- spTransform(bike_racks, CRS=wgs84)
	bike_points <<- SpatialPoints(bike_racks@coords, proj4string=wgs84)

	indices <<- over(bike_points, analysis_area)
	indices[is.na(indices)] = 0
	bike_points = bike_points[(indices == 1),]

	# The first thing to know is that nearly every object in R is really just a list with named elements. 
	# http://www.agapow.net/programming/r/oop-in-r/primitive-r-objects-and-s3/

	bike_pattern = as.ppp(bike_points)
	bike_pattern$window = boundingbox(analysis_area)

	#bike_grid@data[x + (y * width),1]
	bike_density = density(bike_pattern, sigma=0.005, dimyx=c(round(height), width))
	bike_grid = as(bike_density, "SpatialGridDataFrame")
	bike_grid@data = 255.0 * bike_grid@data / max(bike_grid@data)
	writeGDAL(bike_grid, fname="temp.png", drivername="PNG")

	library(rgl)
	bike_matrix = as(bike_grid, "matrix")
	bike_matrix[is.na(bike_matrix)] = 0
	rgl.bg(color = "white")
	rgl.clear()
	rgl.viewpoint(theta = 0, phi = 60, zoom=0.7, scale=c(1,0.2,1))
	rgl.surface(x = 1:nrow(bike_matrix), z = 1:ncol(bike_matrix), y = bike_matrix, 
		texture="/home/michael/mapping-sustainability/analysis_area-osm.png", specular="blue")
	rgl.snapshot("temp3.png")
}


polygon_point_grid = function(polygons, count_colname)
{
	points = NULL

	for (index in 1:(nrow(polygons)))
	{
		point_count = polygons@data[index, count_colname]
		box = bbox(polygons[index,])
		if (!is.na(point_count))
			for (grid in 1:(point_count + 2))
			{
				width = box['x','max'] - box['x','min']
				height = box['y','max'] - box['y','min']

				x = box['x','min'] + (width * ((1:grid) + 1.0) / (grid + 2.0)) 
				y = box['y','min'] + (height * ((1:grid) + 1.0) / (grid + 2.0)) 
				x = rep(x, grid)
				y = unlist(lapply(y, function(x) rep(x,grid)))

				gridpoints = SpatialPoints(data.frame(x,y), proj4string=wgs84)
				indices = over(gridpoints, SpatialPolygons(polygons[index,]@polygons, proj4string=wgs84))
				if (length(indices) > 0)
					indices = indices * (1:length(indices))
				indices = indices[!is.na(indices)]
				# print(paste(index, ": ", grid, "=", length(indices), "of", point_count)) 
				if (length(indices) >= point_count)
				{
					indices = indices[1:point_count]
					if (is.null(points))
						points = gridpoints[indices,]
					else
						points = spRbind(points, gridpoints[indices,])
					break
				}
			}
	}

	if (!is.null(points))
		points = SpatialPointsDataFrame(points@coords, data.frame(id=1:length(points)), proj4string=wgs84)

	return(points)
}					
				
random_barrels = function()
{
	# Create random points to simulate number of rain barrels
	library("plyr")
	library("splancs")

	barrel_zips = readOGR(dsn="/home/michael/mapping-sustainability/census-zip", layer="census-zip")
	barrel_zips <<- spTransform(barrel_zips, CRS=wgs84)

	barrels <<- matrix(ncol=2, nrow=0)
	for (index in 1:nrow(barrel_zips))
		if (!is.na(barrel_zips$RAINBARREL[index]))
		{
			x = csr(barrel_zips@polygons[[index]]@Polygons[[1]]@coords, barrel_zips$RAINBARREL[index])
			barrels <<- rbind(barrels, x)
		}
	row.names(barrels) = 1:nrow(barrels)

	barrels = SpatialPointsDataFrame(coords = barrels, data = data.frame(number = 1:nrow(barrels)), proj4string = wgs84)

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