地理县级热力图

Posted 2023-02-23

【中文标题】地理县级热力图

【英文标题】：Geographical county level heat map

【发布时间】：2017-09-25 09:58:14

【问题描述】：

我经常需要在 R 中创建地理热图。目前，我一直在我的办公室计算机中使用许可版本的 Tableau 来执行此操作，它做得非常好。但是当我不在办公室时，我需要学习如何去做。这些数据有时是机密的，因此我无法通过 Internet 使用 Tableau public。我查看了但找不到任何能产生我需要的结果的解决方案。

数据包含印度贾坎德邦的地区名称以及 6 至 14 岁年龄段的儿童人口（以千计）。在 Tableau 中，我只需将 DISTNAME 列设置为“县”级别的“地理角色”，它就会从互联网 (OpenStreetMap) 中提取州地图以及地区边界，并生成这样的热图，这就是结果如果可能的话，我希望从 R 中获得：

数据是：

geo_data <- structure(list(DISTNAME = c("BOKARO", "CHATRA", "DEOGHAR", "DHANBAD", 
"DUMKA", "GARHWA", "GIRIDIH", "GODDA", "GUMLA", "HAZARIBAGH", 
"JAMTARA", "KHUNTI", "KODARMA", "LATEHAR", "LOHARDAGA", "PAKUR", 
"PALAMU", "PASHCHIMI SINGHBHUM", "PURBI SINGHBHUM", "RAMGARH", 
"RANCHI", "SAHIBGANJ", "SARAIKELA-KHARSAWAN", "SIMDEGA"), POP = c(521.5, 
196.5, 323.8, 445.5, 123, 373.9, 357.6, 248.2, 212.4, 686.7, 
626.7, 383.6, 391.9, 141, 436.1, 454.6, 301.3, 325.5, 193.7, 
238.3, 208.7, 587.4, 130.1, 268)), .Names = c("DISTNAME", "POP"
), row.names = c(NA, 24L), class = "data.frame")

看起来像：

              DISTNAME   POP
1               BOKARO 521.5
2               CHATRA 196.5
3              DEOGHAR 323.8
4              DHANBAD 445.5
5                DUMKA 123.0
6               GARHWA 373.9
7              GIRIDIH 357.6
8                GODDA 248.2
9                GUMLA 212.4
10          HAZARIBAGH 686.7
11             JAMTARA 626.7
12              KHUNTI 383.6
13             KODARMA 391.9
14             LATEHAR 141.0
15           LOHARDAGA 436.1
16               PAKUR 454.6
17              PALAMU 301.3
18 PASHCHIMI SINGHBHUM 325.5
19     PURBI SINGHBHUM 193.7
20             RAMGARH 238.3
21              RANCHI 208.7
22           SAHIBGANJ 587.4
23 SARAIKELA-KHARSAWAN 130.1
24             SIMDEGA 268.0

【问题讨论】：

您需要找到一个 shapefile 才能继续（除非有人知道会自动从 OSM 中提取它们的包）。如果你有一个你感兴趣的县的 shapefile，那么在 R 中映射它们是相当简单的。

您可以查看此问题/答案***.com/questions/7747991/… 以及：computerworld.com/article/3175623/data-analytics/… 正如 Tumbledown 所说，您需要找到适合您所在国家/地区的 shapefile 才能开始使用。

...您应该可以在这里找到这些 shapefile，您只需要将它们过滤到 Jharkhand：gadm.org/download

有字面意思几十个，如果不是数百个例子。你做过任何谷歌搜索吗？如果是这样，什么不起作用？

@hrbrmstr：我进行了搜索，但没有一个产生与我从 Tableau 获得的结果相似的结果。这就是为什么我决定在我的问题中包含 Tableau 结果。

【参考方案1】：

您需要 SHP 文件，可以使用 getData() 找到该文件。完整的工作代码：

library(tidyverse)
library(broom)
library(rgdal)

您的地理数据

geo_data <- structure(list(DISTNAME = c("BOKARO", "CHATRA", "DEOGHAR", "DHANBAD", "DUMKA", "GARHWA", "GIRIDIH", "GODDA", "GUMLA", "HAZARIBAGH", "JAMTARA", "KHUNTI", "KODARMA", "LATEHAR", "LOHARDAGA", "PAKUR", "PALAMU", "PASHCHIMI SINGHBHUM", "PURBI SINGHBHUM", "RAMGARH", "RANCHI", "SAHIBGANJ", "SARAIKELA-KHARSAWAN", "SIMDEGA"),
                       POP = c(521.5, 196.5, 323.8, 445.5, 123, 373.9, 357.6, 248.2, 212.4, 686.7, 626.7, 383.6, 391.9, 141, 436.1, 454.6, 301.3, 325.5, 193.7, 238.3, 208.7, 587.4, 130.1, 268)),
                  .Names = c("DISTNAME", "POP"),
                  row.names = c(NA, 24L),
                  class = "data.frame")

获取地图

library(raster)
IN2 <- getData('GADM', country='IND', level=2)
IN2 <- spTransform(IN2, CRS("+init=epsg:4326"))
IN2_map <- tidy(IN2, region = "NAME_2")

geo_data 中的 id 降低

geo_data$DISTNAME <- tolower(geo_data$DISTNAME)


IN2_map %>%
   mutate(id = tolower(id)) %>%
   left_join(geo_data, by = c("id" = "DISTNAME")) %>%
   ggplot() +
   geom_polygon(aes(long, lat, group=group, fill = POP), color = "black")

【讨论】：

这看起来朝着正确的方向发展。谢谢。现在我必须弄清楚如何更接近使用 Tableau 产生的结果 - 只是带有区边界的州（没有超出州边界），上面写有区名称和填充参数的值。

提示：geom_text

是的，但是geom_text只能在修剪掉州界外的巨大区域后应用。目前，想要的部分显示为图像的一个非常小的区域。应该有一些方法可以放大它。

我发现从国家数据中提取州数据（仅使用基本 R 方法，而不是 data.table 的方法，这对我来说几乎是一种潜意识习惯）然后进行绘图就可以了。跨度>

【参考方案2】：

在下面的解决方案中，我使用了从以下位置下载的地图 shapefile：http://projects.datameet.org/maps/districts/

编辑：后来我还尝试了从http://gadm.org/country 提取的贾坎德邦地图，该地图显示地区边界略有不同。它与互联网上其他州的政治地图更匹配。

这是我的解决方案：

library(tmap)
library(tmaptools)

geo_data <- data.frame(
  DISTNAME = c("BOKARO", "CHATRA", "DEOGHAR", "DHANBAD", "DUMKA", "GARHWA", "GIRIDIH", "GODDA", "GUMLA", "HAZARIBAGH", "JAMTARA", "KHUNTI", "KODARMA", "LATEHAR", "LOHARDAGA", "PAKUR", "PALAMU", "PASHCHIMI SINGHBHUM", "PURBI SINGHBHUM", "RAMGARH", "RANCHI", "SAHIBGANJ", "SARAIKELA-KHARSAWAN", "SIMDEGA"),
  POP = c(521.5, 196.5, 323.8, 445.5, 123, 373.9, 357.6, 248.2, 212.4, 686.7, 626.7, 383.6, 391.9, 141, 436.1, 454.6, 301.3, 325.5, 193.7, 238.3, 208.7, 587.4, 130.1, 268))

# the path to shape file
shp_file <- "H:/Mapping/maps-master/Districts/Census_2011/2011_Dist.shp"

india <- read_shape(shp_file, as.sf = TRUE, stringsAsFactors = FALSE)
india$DISTRICT <- toupper(india$DISTRICT)

jharkhand <- india[india$ST_NM =="Jharkhand", ]

jharkhand_pop <- merge(x = jharkhand,
                       y = geo_data,
                       by.x = "DISTRICT",
                       by.y = "DISTNAME")

#tmap_mode(mode = "plot") # static
tmap_mode(mode = "view") # interactive

qtm(jharkhand_pop, fill = "POP",
    text = "DISTRICT",
    text.size=.9)

静态地图（plot 模式）非常好，但交互式地图（view 模式）非常棒。它提供了从 Internet 的三个不同来源提取其他地图信息的选项。

非常感谢 tmap 和 tmaptools 软件包的创建者。这种方法远胜于网上很多比较冗长和笨拙的解决方案。

如果我们想要更多的定制：

tm_shape(jharkhand_pop) +
  tm_polygons() +
  tm_shape(jharkhand_pop) +
  tm_borders() +
  tm_fill("POP",
          palette = get_brewer_pal("YlOrRd", n = 20),
          n = 20,
          legend.show = F,
          style = "order") + # "cont" or "order" for continuous variable
  tm_text("DISTRICT", size = .7, ymod = .1) +
  tm_shape(jharkhand_pop) +
  tm_text("POP", size = .7, ymod = -.2)

我们在plot 模式下得到以下信息：