Рисуем картограммы в R

Шейп-файлы с картой России (указаны границы административно-территориального деления вплоть до муниципальных районов и городских округов) можно получить на этом сайте сообщества GIS-Lab. Отдельные шейп-файлы для регионов РФ и стран — бывших республик СССР GIS-Lab выкладывает сюда (в файлах boundary-* указаны границы АТО вплоть до районов городских округов и территориальных органов управления).

> # Закажем необходимую для работы библиотеку maptools…
> library("maptools")
> #… и загрузим шейп-файл АТО Удмуртии в объект udm:
> udm <- readShapePoly("~/projects/r-maps/RU-UD/boundary-polygon.shp",
+ proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

Естественно, вместо "~/projects/r-maps/RU-UD/boundary-polygon.shp" следует указать имя Вашего файла, а вместо "+proj=longlat +ellps=WGS84 +datum=WGS84" — Ваши предпочтения по проекции и системе координат (Proj.4 представление). Значение proj4string можно изменять без перезагрузки шейп-файла, например, так:

> proj4string(udm) <- CRS("+proj=merc")

Собственно, уже можно начинать рисовать.

> # Создаём файл-устройство для рисования:
> png("boundary-all.png", width=300, height=400, units="px", bg="transparent")
> # Задаём поля отступов:
> par(mar=c(0.5,0.5,0.5,0.5))
> # Рисуем все АТО в объекте udm:
> plot(udm, col="white")
> # Закрываем файл-устройство:
> dev.off()

Результат:

Естественно, отображением можно управлять. Посмотрим на объект udm подробнее:

> summary(udm)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min      max
x 51.12416 54.43483
y 55.85058 58.54584
Is projected: FALSE 
proj4string :
[+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs]
Data attributes:
     OSM_ID                    NAME       ADMIN_LVL  
 Min.   : -2076064   Октябрьское :  3   8      :311  
 1st Qu.: -1078857   Первомайское:  3   6      : 30  
 Median : -1066632   Каменское   :  2   9      :  5  
 Mean   :   521079   Камское     :  2   10     :  1  
 3rd Qu.: -1017007   Кильмезское :  2   2      :  1  
 Max.   :144692757   (Other)     :341   (Other):  2  
                     NA's        :  4   NA's   :  7

Нам интересны атрибуты OSM_ID, NAME и ADMIN_LVL. Так будет выглядеть та же карта с точностью до муниципальных районов и городских округов:

> png("boundary-6.png", width=300, height=400, units="px", bg="transparent")
> par(mar=c(0.5,0.5,0.5,0.5))
> plot(udm[udm$ADMIN_LVL == 6 & !is.na(udm$ADMIN_LVL), ], col="white")
> dev.off()

Конструкция "!is.na(udm$ADMIN_LVL)" отбрасывает те объекты, у которых не указан административный уровень. Если этого не делать, получим ошибку "NAs not permitted in row index". Запятая перед закрывающей квадратной скобкой (выделена красным) позволяет избежать ошибки "undefined columns selected".

Отбор отображаемых объектов можно делать и по их имени, вот так:

> plot(udm[udm$NAME == "Малопургинский район" & !is.na(udm$NAME), ], col="white")

К сожалению, в данном шейп-файле не содержится сведений о том, какие территориальные органы управления входят в какой район. Эту информацию нужно указывать самостоятельно.