Использовать Rcartogram для объекта SpatialPolygonsDataFrame

Question

Я пытаюсь сделать то же самое, что и в этом вопросе, Картограмма + карта хопоплета в R, но начиная с SpatialPolygonsDataFrame и надеясь в итоге получить объект того же типа.

Я мог бы сохранить объект как шейп- файл, использовать scapetoad, открыть его и преобразовать обратно, но я бы предпочел, чтобы все это было в R, чтобы процедура была полностью воспроизводимой, и чтобы я мог автоматически кодировать десятки вариантов.

Я раздвоил код Rcartogram на github и добавил свои усилия здесь.

По сути, эта демонстрация создает SpatialGrid на карте, ищет плотность населения в каждой точке сетки и преобразует ее в матрицу плотности в формате, необходимом для cartogram() работать на. Все идет нормально.

Но, как интерполировать исходные точки карты на основе выходных данных cartogram()?

Здесь есть две проблемы. Во-первых, чтобы получить карту и сетку в тех же единицах, чтобы позволить интерполяцию. Второе - это получить доступ к каждой точке каждого многоугольника, интерполировать ее и держать все в правильном порядке.

Сетка указана в единицах сетки, а карта в проецируемых единицах (в примере с longlat). Либо сетка должна быть спроецирована в longlat, либо карта в единицах сетки. Моя мысль состоит в том, чтобы сделать поддельный CRS и использовать его вместе с spTransform() функция в package(rgdal), так как это обрабатывает каждую точку объекта с минимальной суетой.

Доступ к каждой точке затруднен, потому что они представляют собой несколько слоев вниз в объект SpPDF: объект> полигоны> полигоны> линии> координаты, я думаю. Любые идеи, как получить к ним доступ при сохранении структуры карты в целом?

Answer 1

Эта проблема может быть решена с помощью getcartr пакет, доступный на GitHub Криса Брансдона, как красиво объяснено в этом сообщении в блоге.

quick.carto Функция делает именно то, что вы хотите - занимает SpatialPolygonsDataFrame в качестве входных данных и имеет SpatialPolygonsDataFrame в качестве вывода.

Воспроизводите суть примера в сообщении в блоге здесь в случае, если ссылка не работает, с моим собственным стилем, смешанным в & исправлены опечатки:

( Shapefile; данные о населении Всемирного банка)

library(getcartr)
library(maptools)
library(data.table)

world <- readShapePoly("TM_WORLD_BORDERS-0.3.shp")
#I use data.table, see blog post if you want a base approach;
#  data.table wonks may be struck by the following step as seeming odd;
#  see here: http://stackru.com/questions/32380338
#  and here: https://github.com/Rdatatable/data.table/issues/1310
#  for some background on what's going on.
world@data <- setDT(world@data)

world.pop <- fread("sp.pop.totl_Indicator_en_csv_v2.csv",
                   select = c("Country Code", "2013"),
                   col.names = c("ISO3", "pop"))

world@data[world.pop, Population := as.numeric(i.pop), on = "ISO3"]

#calling quick.carto has internal calls to the
#  necessary functions from Rcartogram
world.carto <- quick.carto(world, world$Population, blur = 0)

#plotting with a color scale
x <- world@data[!is.na(Population), log10(Population)]
ramp <- colorRampPalette(c("navy", "deepskyblue"))(21L)
xseq <- seq(from = min(x), to = max(x), length.out = 21L)
#annoying to deal with NAs...
cols <- ramp[sapply(x, function(y)
  if (length(z <- which.min(abs(xseq - y)))) z else NA)]

plot(world.carto, col = cols,
     main = paste0("Cartogram of the World's",
                   " Population by Country (2013)"))