R: horoplet map of Russia with an enlarged European part

  • Tutorial

Briefly about the main thing : recently read the post infotanka . I climbed onto Tatyana Misyutina’s website and spied a horoplet map of Russia there with an enlarged European part. And indeed, a really cool idea. Convenient, clear. I wanted to make myself a template for R for the same graphs. After all, good ideas should be propagated by division?

From the website Global Administrative Areas we take the polygons of the borders of the regions of Russia. There just is a native format for R. Data is loaded as a SpatialPolygonsDataFrame .

We load the file, try to draw what is:

rusdf <- load('RUS_adm1.RData’)

It turned out such a strange thing:

In principle, everything is correct. We did not use any projections, in terms of coordinates everything seems to be true. But the "Chukchi scatter" is embarrassing. Of course, the most correct solution would be to use a cylindrical projection to those polygons that already exist. Chukotka will "gather", but the border along the 180 meridian of longitude will remain and will be visible. And on the island of Wrangel, and Chukotka Autonomous Okrug. Apparently, due to errors in one polygon, they do not connect. Poorly. Therefore, it was decided to expand the longitude range over 180 degrees. Yeah, "the number of Pi in wartime can reach four."

for(i in 1:length(gadm@polygons)){
  for(j in 1:length(gadm@polygons[[i]]@Polygons)){
    gadm@polygons[[i]]@Polygons[[j]]@coords[,1]<- sapply(gadm@polygons[[i]]@Polygons[[j]]@coords[,1], function(x){
        if(x < 0){

Something is already looming. But what is it? Chita region, Aginsky Buryat, Koryak and Komi-Permyak autonomous okrugs. How is it in our yard for a year? 2003 or 2013? In general, it is necessary to begin the process of administrative unification of the subjects of the Russian Federation. We will keep the old regions, just in case. Suddenly useful for visualizing archived data.

# Zabaikalsky krai (Chitinskaya obl. + Aginskiy Buryatskiy AOk)
united.reg[united.reg == 2 | united.reg == 13] <- 91
united.reg <- as.character(united.reg)
rus.map <- unionSpatialPolygons(gadm, united.reg)
# Returning old regions (before unioning)
old.regions <- list()
old.regions <- c(old.regions, Polygons(gadm[gadm$ID_1==2,]@polygons[[1]]@Polygons, ID = '2'))
old.regions <- c(old.regions, Polygons(gadm[gadm$ID_1==13,]@polygons[[1]]@Polygons, ID = '13'))
rus.map <- SpatialPolygons(c(slot(rus.map,'polygons'), old.regions)) 

After the unification of the regions, in some places on the site of the old borders, artifacts form due to non-intersection. We filter artifacts according to the small area of ​​the polygons and the hole parameter (function clean.borders). With the preparation of the map, that seems to be all.

You can load a table with statistics and draw.

The table shows the IDs of regions for which polygons will be combined with statistics. There are outdated regions, and existing now. The polygons for rendering are selected according to the data column: if “NA” is in the statistics field, this subject of the Federation is not rendered at all; if the statistics are “zero”, the region is grayed out. The remaining colors of the map are configured using the RColorBrewer palettes. Here you can choose the type of palette: (1) sequential - shades of one color showing the severity of the sign or (2) diverging - deviation from the average value in plus (one color, e.g. green) or minus (another color, e.g. red).

The next step is to create a basemap. The subsequent summary picture will simply be a combination of the two views of this base chart. We remove from the base chart all indents, background, axis labels and labels.

p <- p + theme(axis.line=element_blank(),axis.text.x=element_blank(),
               legend.position = 'none',
               panel.margin = unit(c(0,0,0,0), 'cm'),
               axis.ticks.margin = unit(0, 'cm'),
               axis.ticks.length = unit(0.001, 'cm'),
               plot.margin = unit(c(0,0,0,0), 'cm'),
               panel.grid = element_blank(),
               panel.background = element_blank()
p <- p + labs(x=NULL, y = NULL)

And a small “hack” related to the order of drawing polygons. When building polygons from SpatialPolygons, Ggplot2 does not take into account their hole property. As a result, if the Moscow region will be drawn after Moscow, even though the source data contained an “exclusive” hole polygon, ggplot will fill the territory of Moscow with the color of the region. The same problem with Adygea. Therefore, these two regions are redrawn last, after all the others are built.

Now from one graph we need to get two views: the enlarged European part (p1) and the rest of the country (p2). At the same stage, select the map projection (with equal distances), the position of the "camera" and the viewing angle. We return the legend to one of the graphs. The beauty of the coord_map function is that when rendering data that does not appear on the picture is not filtered out and taken into account when constructing the visible part of the graph. That is, even if the range of statistics for the European and Asian parts of Russia is different, there will be no glitches in the legend and color coding. I selected the fields and viewing angles for this graph using the “scientific poking method”, sorting through various options. I would be glad if someone tells me a more reliable and repeatable method.

The final picture is assembled using the grid library. It allows you to mark out areas for inserting individual elements into the overall picture (viewport's). So far, without “heavy” graphics, with the help of placeholders we

’ll figure out a suitable option for placing elements: We draw a simple icon of the magnifying glass to show the difference in the scale of the left and right parts of the graph.

magnif.glass <- function(vport){
  grid.circle(x=.6,y=.6,r=.3, gp=gpar(lwd=1.5, col='grey70'), vp = vport)
  grid.lines(x=c(.6,.6), y=c(.5,.7), gp=gpar(lwd=1.5, col='grey70'), vp = vport)
  grid.lines(x=c(.5,.7), y=c(.6,.6), gp=gpar(lwd=1.5, col='grey70'), vp = vport)
  grid.lines(x=c(.1,.4), y=c(.1,.4), gp=gpar(lwd=1.5, col='grey70'), vp = vport)
  grid.lines(x=c(.1,.3), y=c(.1,.3), gp=gpar(lwd=3, col='grey70'), vp = vport)

And we add a text leader with an average value or some other huge and important figure. That seems to be all.

All code with test data lies on GitHub . At the end, two questions to the community:
1. Where can I get a shape-file of the borders of the new Moscow? I would like to update this moment.
2. What exactly needs to be redone / improved? And then I'm a fake welder. Tell me?

Also popular now: