library("rgdal")

# Convert LEED locations CSV to a shapefile

csv = read.csv("leed-locations.csv", colClasses = "character", fill=T)
print(paste(nrow(csv), "total locations"))

# Create points from valid long/lat

csv[,"latitude"] = as.numeric(csv[,"latitude"])
csv[,"longitude"] = as.numeric(csv[,"longitude"])
csv = csv[(!is.na(csv[,"latitude"]) & !is.na(csv[,"latitude"]) & 
	(csv[,"latitude"] < 180) & (csv[,"latitude"] > -180) &
	(csv[,"latitude"] != 0) & (csv[,"latitude"] != 0)),]
csv_points = SpatialPoints(coords = csv[, c("longitude","latitude")], 
	proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

print(paste(nrow(csv), "filtered locations"))

# Find analysis cities

analysis_areas <<- readOGR(dsn="/home/michael/mapping-sustainability/boundaries/analysis-areas",
	layer="boundaries-county")
analysis_areas <<- spTransform(analysis_areas, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
analysis_polygons <<- SpatialPolygons(analysis_areas@polygons, proj4string=CRS(proj4string(analysis_areas)))
indices = over(csv_points, analysis_polygons)
csv$City = as.character(analysis_areas@data[indices, "CITY"])


# Create data frame containing only points within analysis cities

geospatial = SpatialPointsDataFrame(coords = csv_points, data = csv)
geospatial = geospatial[!is.na(geospatial@data$City),]


# Save as shapefile

writeOGR(geospatial, dsn="/home/michael/mapping-sustainability/green-buildings",
	layer="leed-locations", driver="ESRI Shapefile", overwrite_layer=T)
