r - 在 SpatialPolygonsDataFrame 对象上使用 Rcartogram

Question

我正在尝试做与这个问题相同的事情，Cartogram + choropleth map in R ，但从 SpatialPolygonsDataFrame 开始并希望以相同类型的对象结束。

我可以将对象保存为 shapefile，使用 scapetoad ，重新打开它并转换回来，但我宁愿将它全部放在 R 中，以便该过程完全可重现，并且我可以自动编码数十种变体。

我已经在 github 上 fork 了 Rcartogram 代码并添加了我迄今为止的努力 here .

本质上，这个演示所做的是在 map 上创建一个 SpatialGrid，查找网格每个点的人口密度并将其转换为 cartogram() 所需格式的密度矩阵。工作。到现在为止还挺好。

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

这里有两个问题。第一种是将 map 和网格放入相同的单位以允许插值。第二个是访问每个多边形的每个点，对其进行插值，并保持它们的顺序正确。

网格以网格单位为单位， map 以投影单位为单位(在示例 longlat 的情况下)。要么必须将网格投影到 longlat 中，要么将 map 投影到网格单位中。我的想法是制作一个假的 CRS 并将其与 spTransform() 一起使用。函数在 package(rgdal) ，因为这会以最小的麻烦处理对象中的每个点。

访问每个点都很困难，因为它们在 SpPDF 对象中有几层:对象>多边形>多边形>线条>我认为坐标。任何想法如何在保持整体 map 结构完整的同时访问这些？

Answer 1

这个问题可以用 getcartr 解决包装，可在 Chris Brunsdon's GitHub 上获得，正如 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)"))