7

我正在尝试做与此问题中提出的相同的事情，即 Cartogram + choropleth map in R，但从 SpatialPolygonsDataFrame 开始并希望最终得到相同类型的对象。

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

我已经在 github 上分叉了 Rcartogram 代码，并在此处添加了我的努力。

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

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

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

网格采用网格单位，地图采用投影单位（在示例 longlat 的情况下）。必须将网格投影到 longlat，或者将地图投影到网格单元中。我的想法是制作一个假的 CRS 并将其与 in 中的spTransform()函数一起使用package(rgdal)，因为它可以轻松处理对象中的每个点。

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

4

1 回答 1

2

这个问题可以通过Chris Brunsdon 的 GitHubgetcartr上的包来解决，正如这篇博文中所阐述的那样。

该quick.carto函数完全符合您的要求——将 aSpatialPolygonsDataFrame作为输入并以 aSpatialPolygonsDataFrame作为输出。

如果链接失效，请在此处复制博客文章中示例的精髓，并混合我自己的风格并修复错别字：

（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)"))

在此处输入图像描述

于 2015-09-08T23:43:52.260 回答