在SpatialPolygonsDataFrame对象上使用Rcartogram

Question

我正在尝试做这个问题中问到的相同的事情，Cartogram + choropleth map in R，但是从一个SpatialPolygonsDataFrame开始，并希望得到相同类型的对象。

我可以将对象保存为shapefile，使用scapetoad，重新打开它并转换回来，但我更愿意将其全部放在R中，这样过程就可以完全重现，这样我就可以自动编写几十种变体。

我已经在github上派生了Rcartogram代码，并添加了到目前为止here的工作。

本质上，此演示所做的是在地图上创建SpatialGrid，查找网格中每个点的人口密度，并将其转换为cartogram()所需格式的密度矩阵。到目前一切尚好。

但是，如何根据cartogram()的输出对原始地图点进行插值？

这里有两个问题。第一种方法是将地图和网格放在相同的单元中，以便进行插值。第二种方法是访问每个多边形的每个点，对其进行插值，并使它们保持正确的顺序。

栅格使用栅格单位，而贴图使用投影单位(在示例longlat的情况下)。要么将网格投影到longlat中，要么将地图投影到网格单元中。我的想法是创建一个假的CRS，并将其与package(rgdal)中的spTransform()函数一起使用，因为这样可以最大限度地处理对象中的每个点。

访问每个点都很困难，因为它们位于SpPDF对象的几层之下:我认为是object>polygons>Polygons>lines>coords。你知道如何在保持整体地图结构完整的情况下访问这些地图吗？

Answer 1

这个问题可以通过Chris Brunsdon's GitHub上的getcartr包来解决，this博客文章中对此进行了很好的解释。

quick.carto函数做的正是您想要的--接受SpatialPolygonsDataFrame作为输入，并将SpatialPolygonsDataFrame作为输出。

在这里重现博客文章中示例的精髓，以防链接死掉，我自己的风格混合在一起&拼写错误修复了：

(Shapefile；World Bank population data)

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://stackoverflow.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)"))

​

​