Plot Chicago ward map with labels that are actually in the wards

wardmap.R

##==============================================================================
## 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")))