Skip to content

Instantly share code, notes, and snippets.

@geneorama

geneorama/wardmap.R

Last active Jan 2, 2020
Embed
What would you like to do?
Plot Chicago ward map with labels that are actually in the wards
##==============================================================================
## INITIALIZE
##==============================================================================
rm(list=ls())
library(geneorama) ## Not actually needed
library(data.table)
library("rgeos")
library("leaflet")
library("colorspace")
library("sp")
# library("spdep")
# library("rgdal")
# library("RColorBrewer")
# library("ggplot2")
# library("bit64")
## Source for map coloring algo - https://github.com/hunzikp/MapColoring
# devtools::install_github("hunzikp/MapColoring")
library(MapColoring)
## Define globals
# LEAFLET_TILES <- "Stamen.TonerHybrid"
LEAFLET_TILES <- "OpenStreetMap.Mapnik"
##==============================================================================
## DOWNLOAD DATA
##==============================================================================
## Contained code to download Chicago's wards
shpCityWards <- local({
cur <- getwd()
on.exit(setwd(cur))
tmp <- tempfile(fileext = ".zip")
setwd(dirname(tmp))
# url <- "https://data.cityofchicago.org/api/geospatial/sp34-6z76?method=export&format=Shapefile"
url <- "https://data.cityofchicago.org/api/geospatial/sp34-6z76?method=export&format=GeoJSON"
download.file(url, destfile = tmp)
shp <- rgdal::readOGR(basename(tmp), stringsAsFactors = FALSE)
shp
})
## Generate city outline
shpCityOutline <- rgeos::gUnaryUnion(as(shpCityWards, "SpatialPolygons"))
##==============================================================================
## CALCULATE WARD COLORS
##==============================================================================
## Get colors for each ward using 4 color algorithm
pal <- colorspace::qualitative_hcl(n = 4, h = c(26, -264), c = 70, l = 70)
colors <- pal[MapColoring::getColoring(shpCityWards)]
## Plot to check if colors worked
plot(shpCityWards, col=colors)
##==============================================================================
## LABEL WITH DEFAULT LABEL LOCATIONS
##==============================================================================
## Extract labels by pulling "labpt"
## This is how you *should* be able to do it
# labs <- sapply(shpCityWards@polygons, `@`, "labpt")
## Extract labels by pulling "labpt"
labs <- vector("list", length(shpCityWards@polygons))
for(i in 1:length(shpCityWards@polygons)){
labs[[i]] <- c(shpCityWards@polygons[[i]]@labpt,
as.numeric(shpCityWards@polygons[[i]]@ID))
}
labs <- data.table(do.call(rbind, labs))
setnames(labs, c("x","y", "ward"))
## Plot in leaflet
leaflet() %>%
addProviderTiles(LEAFLET_TILES) %>%
addPolygons(data = shpCityOutline, fill = FALSE, color = "black", weight = 5) %>%
addPolygons(data = shpCityWards, fillColor = colors, fillOpacity = 0.5,
weight = 0.5, label = ~paste("Ward:", ward) ) %>%
addLabelOnlyMarkers(data = labs, ~x, ~y, label = ~as.character(ward),
labelOptions = labelOptions(noHide = TRUE,
direction = "center",
offset = c(0, 0), opacity = 1,
textsize = "12px", textOnly = TRUE,
style = list("font-style" = "bold")))
##==============================================================================
## LABEL WITH gCentroid LOCATIONS
##==============================================================================
## That wasn't too great.
## Try with creating labs with gCentroid
## Create labels
labs <- as.data.frame(gCentroid(shpCityWards, byid = TRUE))
labs$ward <- shpCityWards$ward
## Plot in leaflet
leaflet() %>%
addProviderTiles(LEAFLET_TILES) %>%
addPolygons(data = shpCityOutline, fill = FALSE, color = "black", weight = 5) %>%
addPolygons(data = shpCityWards, fillColor = colors, fillOpacity = 0.5,
weight = 0.5, label = ~paste("Ward:", ward) ) %>%
addLabelOnlyMarkers(data = labs, ~x, ~y, label = ~as.character(ward),
labelOptions = labelOptions(noHide = TRUE,
direction = "center",
offset = c(0, 0), opacity = 1,
textsize = "12px", textOnly = TRUE,
style = list("font-style" = "bold")))
## That was better, but still way off for wards like 2, 27, 29, 41, etc
##==============================================================================
## LABEL WITH centroid FUNCTION
##==============================================================================
## Source for functions:
## https://gis.stackexchange.com/a/265475/78424
#' find the center of mass / furthest away from any boundary
#'
#' Takes as input a spatial polygon
#' @param pol One or more polygons as input
#' @param ultimate optional Boolean, TRUE = find polygon furthest away from centroid. False = ordinary centroid
require(rgeos)
require(sp)
centroid <- function(pol,ultimate=TRUE,iterations=5,initial_width_step=10){
if (ultimate){
new_pol <- pol
# For every polygon do this:
for (i in 1:length(pol)){
width <- -initial_width_step
area <- gArea(pol[i,])
centr <- pol[i,]
wasNull <- FALSE
for (j in 1:iterations){
if (!wasNull){ # stop when buffer polygon was alread too small
centr_new <- gBuffer(centr,width=width)
# if the buffer has a negative size:
substract_width <- width/20
while (is.null(centr_new)){ #gradually decrease the buffer size until it has positive area
width <- width-substract_width
centr_new <- gBuffer(centr,width=width)
wasNull <- TRUE
}
# if (!(is.null(centr_new))){
# plot(centr_new,add=T)
# }
new_area <- gArea(centr_new)
#linear regression:
slope <- (new_area-area)/width
#aiming at quarter of the area for the new polygon
width <- (area/4-area)/slope
#preparing for next step:
area <- new_area
centr<- centr_new
}
}
#take the biggest polygon in case of multiple polygons:
d <- disaggregate(centr)
if (length(d)>1){
biggest_area <- gArea(d[1,])
which_pol <- 1
for (k in 2:length(d)){
if (gArea(d[k,]) > biggest_area){
biggest_area <- gArea(d[k,])
which_pol <- k
}
}
centr <- d[which_pol,]
}
#add to class polygons:
new_pol@polygons[[i]] <- remove.holes(new_pol@polygons[[i]])
new_pol@polygons[[i]]@Polygons[[1]]@coords <- centr@polygons[[1]]@Polygons[[1]]@coords
}
centroids <- gCentroid(new_pol,byid=TRUE)
}else{
centroids <- gCentroid(pol,byid=TRUE)
}
return(centroids)
}
#Given an object of class Polygons, returns
#a similar object with no holes
remove.holes <- function(Poly){
# remove holes
is.hole <- lapply(Poly@Polygons,function(P)P@hole)
is.hole <- unlist(is.hole)
polys <- Poly@Polygons[!is.hole]
Poly <- Polygons(polys,ID=Poly@ID)
# remove 'islands'
max_area <- largest_area(Poly)
is.sub <- lapply(Poly@Polygons,function(P)P@area<max_area)
is.sub <- unlist(is.sub)
polys <- Poly@Polygons[!is.sub]
Poly <- Polygons(polys,ID=Poly@ID)
Poly
}
largest_area <- function(Poly){
total_polygons <- length(Poly@Polygons)
max_area <- 0
for (i in 1:total_polygons){
max_area <- max(max_area,Poly@Polygons[[i]]@area)
}
max_area
}
labs <- centroid(pol = shpCityWards,
ultimate = TRUE,
iterations = 150,
initial_width_step = .01)
warnings()
labs$ward <- shpCityWards$ward
leaflet() %>%
addProviderTiles(LEAFLET_TILES) %>%
addPolygons(data = shpCityOutline, fill = FALSE, color = "black", weight = 5) %>%
addPolygons(data = shpCityWards, fillColor = colors, fillOpacity = 0.5,
weight = 0.5, label = ~paste("Ward:", ward) ) %>%
addLabelOnlyMarkers(data = labs, ~labs$x, ~labs$y, label = ~as.character(ward),
labelOptions = labelOptions(noHide = TRUE,
direction = "center",
offset = c(0, 0), opacity = 1,
textsize = "12px", textOnly = TRUE,
style = list("font-style" = "bold")))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.