在一个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。有没有想法如何访问这些并保持整个地图结构的完整性？

- duncandoo

---

我在发布了我的问题并且自己使用Rcartogram时，偶然发现了这个问题。到目前为止，我的建议是使用ScapeToad；我正在尝试决定是否可以将其简单性移植到R中。 - MichaelChirico

Answer 1

这个问题可以通过使用可在Chris Brunsdon的GitHub上获得的getcartr包来解决，如此博客文章所详细阐述。 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://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)"))

- MichaelChirico