R: choir map of Russia with a larger European part

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



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



We load the file, we try to draw what we have:

rusdf <- load('RUS_adm1.RData') plot(gadm) 

It turned out such a strange thing:



')

In principle, everything is correct. We didn’t use any projections; in coordinates, everything seems to be true. But the "Chukchi spread" confuses. Of course, the most correct solution would be to use a cylindrical projection to those polygons that already exist. Chukotka will “gather,” but the longitude along the longitude 180 will remain and will be visible. And on Wrangel Island, and on the Chukotka Autonomous Region. Apparently, because of the errors in one polygon, they are not connected. Poorly. Therefore, it was decided to expand the range of longitudes over 180 degrees. Yeah, "Pi may be up to four in wartime."

 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){ x<-359.999+x } else{x} }) } } 

Something is already looming. But what is it? The Chita region, the Aginsk Buryat, Koryak and Komi-Perm autonomous districts. How is our year in the yard? 2003 or 2013? In general, it is necessary to begin the process of administrative association of the subjects of the Russian Federation. The old regions, too, just in case, save. 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, artifacts are formed here and there at the place of the old borders because of non-intersection. We filter artifacts by a small area of ​​the polygons and the hole parameter (the clean.borders function). With the preparation of the card, it seems everything.



You can load a table with statistics and draw.



The table shows the IDs of the regions for which the polygons will be combined with statistics. There are both obsolete regions and existing ones. The selection of polygons for drawing takes place in the data column: if “NA” is in the statistics field, this Federation subject is not drawn at all; if the statistics is “zero” - the region is drawn in gray. The rest of the map colors are configured using RColorBrewer palettes. Here you can also select the type of palette: (1) sequential - shades of the same color, indicating the severity of the sign or (2) diverging - deviation from the mean value by plus (one color, eg green) or minus (another color, eg red).



The next step is to create a base map. The subsequent final picture will be just a combination of two types on this baseline. We remove all indents, backgrounds, axis labels and labels from the base graph.

 p <- p + theme(axis.line=element_blank(),axis.text.x=element_blank(), axis.text.y=element_blank(),axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=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" associated with the order of drawing polygons. Ggplot2 when constructing polygons from SpatialPolygons does not take into account their hole property. As a result, if the Moscow region will be drawn after Moscow, even though there was an "exclusion" polygon-hole in the source data - ggplot will paint over 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: an enlarged European part (p1) and the rest of the country (p2). At the same stage, we select the map projection (with equal distances), the position of the “camera” and the angle of view. On one of the graphs return the legend. The beauty of the coord_map function is that when drawing data that is not on the picture, it is not filtered out and is taken into account when building the visible part of the graph. That is, even if the range of statistics values ​​for the European and Asian parts of Russia will differ - there will be no failures in the legend and color coding. I selected the fields and viewing angles for this chart using the “method of scientific typing”, going through various options. I would be glad if someone tells you a more reliable and repeatable way.



Assembling the final image is done using the grid library. It allows you to mark the area for insertion of individual elements in the overall picture (viewports). So far, without “heavy” graphics, with the help of place holders we pretend a suitable placement option for the elements:







We draw a simple magnifying glass icon 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 add a text callout with an average value or some other huge and important number. That seems to be all.







All code with test data lies on GitHub . At the end there are two questions to the community:

1. Where to get the shape-file of the borders of the new Moscow? I want to update this moment.

2. What exactly needs to be redone / improved? Because I’m a fake welder. Tell me?