Картограммы в R: меняем проекцию

Для смены проекции карты в R можно использовать библиотеку rgdal:

> install.packages("rgdal")
> library("rgdal")

В Ubuntu для работы этой библиотеки требуются пакеты libgdal1-dev и libproj-dev с их зависимостями.

Мы можем скачать и прочитать шейп-файл, как указано ранее (это, кстати, существенно быстрее), однако продемонстрирую вариант подключения к открытой базе геоданных PostGIS при помощи того же rgdal:

> # Читаем список доступных слоёв базы GEN…
> ogrListLayers(dsn="PG:host=gis-lab.info user=guest dbname=gen password=guest")
 [1] "modis_tiles"              "oblasts"                  "ru_adm2_country"
 …
[26] "wrs2_landsat_rus"
> # … и загружаем слой границ России в объект rus.
> rus <- readOGR(dsn="PG:host=gis-lab.info user=guest dbname=gen password=guest",
+ layer="ru_adm2_country", verbose=T)

Проект GIS-Lab для отображения Сибирской части России и для отображения России целиком пользуется равновеликой конической проекцией Альберса с главными параллелями 52° и 64° с.ш. и главным меридианом 105° в.д.. Проекция объекта rus —

> proj4string(rus)
[1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

— что делает его изображение разорванным по 180-му меридиану, в районе Чукотки. Сменим проекцию и нарисуем карту:

> rus2 <- spTransform(rus, CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105
+ +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs"))
> plot(rus2, col="#31A35488", border=F)

Чтобы наложить на полученную карту сетку координат, необходимо эту сетку специально подготовить:

> grd <- gridlines(rus,
+ norths = seq(10, 90, by=5), # рисуем линии каждые 5° между 10 и 90° с.ш.
+ # Отобразятся только те из них, которые попадут в область рисуемой карты.
+ easts = seq(-179.9999, 180, by=5), # Если указать -180, этот меридиан не отобразится.
+ ndiscr=360) # Число опорных точек — если указать мало,
> # координатная сетка после перепроецирования может получиться ломаной.
> # Меняем проекцию координатной сетки так же, как у карты…
> grd2 <- spTransform(grd, CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105
+ +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs"))
> # … и добавляем её к рисунку.
> plot(grd2, add=T, col="#08519C88", lty=2)

Результат приведен в начале заметки.

Перед сменой проекции мы можем упростить контур объекта rus,

> rus.thiny <- thinnedSpatialPoly(rus,
+ tolerance=0.5, # Интенсивность упрощения, чем больше, тем проще контур.
+ minarea=0.1) # Минимальная площадь сохранямых островов.

объем данных уменьшается очень существенно, но машинное время, затраченное на этот процесс, превышает время на смену проекции и отрисовку исходного объекта вместе взятое. Результат с упрощенным контуром — ниже.