Tilegram of Russian regions

README.md

This is a tilegram of Russian regions with area being distorted according to 2015 population. I relied on Gastner and Newman, 2004 distortion algorithm. Each hexagon is approximately equal to 25,000 people. Hover over the regions to see their names and population counts (in Russian).

# (c) Dmitriy Skougarevskiy, November 2016
#
# Use of this source code is governed by the MIT license
# located at https://opensource.org/licenses/MIT

# Load dependencies
packages <- c("maptools", "raster", "sp", "rgeos", "devtools")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
	install.packages(setdiff(packages, rownames(installed.packages())))  
}
lapply(packages, library, character.only = TRUE)
# NB: uncomment the below two packages lines to install
# the packages from Github. They require http://www.fftw.org
# (e.g. brew install fftw)
#devtools::install_github("omegahat/Rcartogram")  
#devtools::install_github("chrisbrunsdon/getcartr", subdir = "getcartr")
library(Rcartogram)
library(getcartr)

# Function to download and open zipped shapefile
# via http://thebiobucket.blogspot.ch/2013/09/batch-downloading-zipped-shapefiles.html
url_shp_to_spdf <- function(URL) {

	require(maptools)

	wd <- getwd()
	td <- tempdir()
	setwd(td)

	temp <- tempfile(fileext = ".zip")
	download.file(URL, temp)
	unzip(temp)

	shp <- dir(tempdir(), "*.shp$")
	y <- readShapePoly(shp)

	unlink(dir(td))
	setwd(wd)
	return(y)
}

# Download Russian regional boundaries shapefile
# that has 2015 regional population embedded therein.
# I aggressively simplified it with mapshaper and by hand
# (removal of the Far North islands, all small islands or lakes
# not exceeding than Sebastopol in area)
rus_bounds <- url_shp_to_spdf("http://donutholes.ch/test/rus_regions_simpl.shp.zip")

# Use Albers Siberia projection
WGS_84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
Albers_Siberia <- CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105 +x_0=18500000 +y_0=0 +ellps=krass +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12 +units=m +no_defs")

proj4string(rus_bounds) <- WGS_84
rus_bounds <- spTransform(rus_bounds, Albers_Siberia)

# Create the cartogram (takes time and CPU, decrease the resolution
# to 512 to enjoy reasonable computation time at the cost of larger error in shape areas)
rus_bounds_transf <- quick.carto(rus_bounds, rus_bounds$region_pop, res = 3000, gapdens = 2, blur = 0)

# A hack to fix possible topology errors
# (via http://gis.stackexchange.com/a/163480)
rus_bounds_transf <- gBuffer(rus_bounds_transf, byid = TRUE, width = -0.0001)

# Export transformed shape now if you feel ready to explore it
#writePolyShape(rus_bounds_transf, "rus_bounds_transf")

proj4string(rus_bounds_transf) <- Albers_Siberia

# Check whether areas are now proportional to population
## Correlation between region population and area now exceeds 0.9:
## cor(rus_bounds_transf@data$region_pop, gArea(rus_bounds_transf, byid = T))
## Difference in area between Crimea and Moscow now almost corresponds to difference in population:
## gArea(rus_bounds_transf, byid = T)[36]/gArea(rus_bounds_transf, byid = T)[54]
## rus_bounds_transf@data$region_pop[36]/rus_bounds_transf@data$region_pop[54]

# Also examine the residuals of the transformation
## discrep <- cart.resids(rus_bounds_transf, rus_bounds$region_pop)
## shs <- auto.shading(c(discrep,-discrep), n = 5, col = brewer.pal(5, "RdBu"))
## choropleth(rus_bounds_transf, discrep, shading = shs)
## choro.legend(600, 600, shs)

# Establish the average number of people per unit of transformed shapefile area
pop_per_area <- mean(rus_bounds_transf@data$region_pop/gArea(rus_bounds_transf, byid = T))
# Set area for hexagons so that each one encompasses 25000 people
hex_area <- 25000/pop_per_area

# Function to produce hexagonal grid
# (via http://strimas.com/spatial/hexagonal-grids/)
make_grid <- function(x, cell_diameter, cell_area, clip = FALSE) {
	if (missing(cell_diameter)) {
		if (missing(cell_area)) {
			stop("Must provide cell_diameter or cell_area")
		} else {
			cell_diameter <- sqrt(2 * cell_area / sqrt(3))
		}
 	}
	ext <- as(extent(x) + cell_diameter, "SpatialPolygons")
	projection(ext) <- projection(x)
	# generate array of hexagon centers
	g <- spsample(ext, type = "hexagonal", cellsize = cell_diameter, 
					offset = c(0.5, 0.5))
	# convert center points to hexagons
	g <- HexPoints2SpatialPolygons(g, dx = cell_diameter)
	# clip to boundary of study area
	if (clip) {
		g <- gIntersection(g, x, byid = T)
	} else {
		g <- g[x, ]
	}
	# clean up feature IDs
	row.names(g) <- as.character(1:length(g))
	return(g)
}

# Create the hexagonal grid
source_hex_grid <- make_grid(rus_bounds_transf, cell_area = hex_area, clip = F)

# Add attributes to the grid to produce the final object
rus_hex <- SpatialPolygonsDataFrame(source_hex_grid, over(source_hex_grid, rus_bounds_transf))

proj4string(rus_hex) <- Albers_Siberia

# Define the color palette
region_colors <- data.frame(REGION = levels(rus_hex$REGION),
							color = brewer.pal(12, name = "Set3")[sample(1:12, 85, replace = T)]
							)

# Add color palette to the grid
rus_hex$color <- as.character(region_colors[match(rus_hex$REGION, region_colors$REGION), "color"])

# Export the resulting shapefile with hexagons
writePolyShape(rus_hex, "rus_hex")

# Convert shapefile to TopoJSON with
# https://github.com/topojson/topojson-1.x-api-reference/blob/master/Command-Line-Reference.md
system(paste0("cd ", getwd(), "; ",
		"topojson --width 960 --height 550 --margin 10 -s .25 --shapefile-encoding utf8 --id-property SP_ID -p oktmo=OKTMO -p region=REGION -p population=region_pop -p color=color -o rus_hex.topojson rus_hex.shp"
		)
	)

index.html

<!DOCTYPE html>
<meta charset="utf-8">
<style>
	.region {
		fill: #fff;
		stroke: #fff;
		stroke-width: 0.1px;
	}
	.region_boundary {
		fill: none;
		stroke: #fff;
		stroke-width: 1px;
	}
	.region:hover {
		fill: #666;
	}
	.hidden {
		display: none;
	}
	div.tooltip {
		color: #222;
		padding: .5em;
		font-family: "Helvetica", Arial, Serif;
		font-size: 10px;
		text-shadow: #f5f5f5 0 1px 0;
		position: absolute;
	}
</style>

<body>
	<script src="http://d3js.org/d3.v3.min.js"></script>
	<script src="http://d3js.org/topojson.v1.min.js"></script>
	<script>
		var width = 960,
			height = 550;

		// We use a pre-projected JSON
		var path = d3.geo.path()
			.projection(null);
	
		var svg = d3.select("body").append("svg")
			.attr("width", width)
			.attr("height", height)
			.append("g");
			
		var tooltip = d3.select('body').append('div')
			.attr('class', 'hidden tooltip');
			
		d3.json('http://donutholes.ch/test/rus_hex.topojson', function(err, input) {
			if (err) return console.error(err);
			
			// Read the TopoJSON structure
			var rus_hex = topojson.feature(input, input.objects.rus_hex).features;

			var g = svg.selectAll('.region')
				.data(rus_hex).enter();
			
			g.append('path')
				.style('fill', function(d) { return d.properties.color; })
				.attr('class', 'region')

				.attr('d', path)
				.on('mousemove', function(d) {
					var mouse = d3.mouse(svg.node()).map(function(d) {
						return parseInt(d);
					});
					tooltip.classed('hidden', false)
						.attr('style', 'left:' + (mouse[0] + 15) +
								'px; top:' + (mouse[1] - 35) + 'px')
						.html(d.properties.region + ",<br> население: " + d3.format(".3s")(d.properties.population));
				})
				.on('mouseout', function() {
					tooltip.classed('hidden', true);
				});
			
			svg.append('path')
    			.datum(topojson.mesh(input, input.objects.rus_hex, function(a, b) { return a.properties.oktmo !== b.properties.oktmo; }))
    			.attr('class', 'region_boundary')
    			.attr('d', path);
						
		});


</script>
</body>