用ggplot制作Cloropleth

huangapple
go评论94阅读模式
英文:

Cloropleth using ggplot

问题

我一直在尝试绘制一个色彩地图。我取得了一些进展，但现在卡住了。这个想法是在巴西地图上标记每个州的暴力故意死亡数量。

我正在使用以下数据集：

  1. State <- c("SP", "SC", "DF", "MG", "RS", "MS", "PR", "AC", "PI", "TO", "MT", "RO", "GO", "RJ", "ES", "MA", "PB", "AL", "RN", "PA", "SE", "PE", "RR", "CE", "AM", "BA", "AP")
  2. VID <- c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8)
  3. VIDBrazil <- data.frame(State, VID)

我从这里获取巴西的geoJSON数据文件。

我正在尝试使用的代码基于这里找到的资源，特别是这里，这里，这里和这里。

我正在使用以下代码：

  1. State <- c("SP", "SC", "DF", "MG", "RS", "MS", "PR", "AC", "PI", "TO", "MT", "RO", "GO", "RJ", "ES", "MA", "PB", "AL", "RN", "PA", "SE", "PE", "RR", "CE", "AM", "BA", "AP")
  2. VID <- c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8)
  3. VIDBrazil <- data.frame(State, VID)
  5. spdf_Brasil <- geojson_read("C:/Users/Giovanni/Desktop/R/brazil_geo.json", what = "sp") #导入JSON文件
  7. spdf_Brasil_A <- tidy(spdf_Brasil) #转换为数据框
  9. ggplot() + #绘制地图
  10.   geom_polygon(data = spdf_Brasil_A, aes( x = long, y = lat, group = group), fill="#69b3a2", color="white") +
  11.   theme_void() +
  12.   coord_map()
  14. spdf_Brasil_MVI <- spdf_Brasil_A %>% #添加用于地图颜色编码的死亡值
  15.   left_join(. , VIDBrazil, by=c("id"="State"))
  17. ggplot() + #带颜色编码的绘图
  18.   geom_polygon(data = spdf_Brasil_MVI, aes( x = long, y = lat, group = group, fill= VID), color="green") + 
  19.   theme_void() +
  20.   scale_fill_gradient(name = "Mortes por 100 mil habitantes", low = "#fdc0b2", high = "#9a1f03", limits = c(0, 60.0)) +
  21.   coord_map()

获得的地图

我缺少什么？我应该怎么做才能得到地图的颜色编码？

谢谢你的帮助！

英文:

I've been trying to plot a cloropleth map. I've made some progress, but I'm stuck now.
The idea is to plot the map of Brazil and color code each state bases on the number of violent intentional deaths in 2021.

I'm using the following dataset:

  1. State &lt;- c(&quot;SP&quot;, &quot;SC&quot;, &quot;DF&quot;, &quot;MG&quot;, &quot;RS&quot;, &quot;MS&quot;, &quot;PR&quot;, &quot;AC&quot;, &quot;PI&quot;,&quot;TO&quot;, &quot;MT&quot;, &quot;RO&quot;, &quot;GO&quot;, &quot;RJ&quot;, &quot;ES&quot;, &quot;MA&quot;, &quot;PB&quot;, &quot;AL&quot;, &quot;RN&quot;, &quot;PA&quot;, &quot;SE&quot;, &quot;PE&quot;, &quot;RR&quot;, &quot;CE&quot;, &quot;AM&quot;, &quot;BA&quot;, &quot;AP&quot;)
  2. VID &lt;- c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8)
  3. VIDBrazil &lt;- data.frame(State, VID)

The geoJSON data file from Brazil I'm sourcing from this Brasil JSON file

The code I'm trying to use is based on resources found here, especially here, here, here and here

The code I'm using is the following:

  1. State &lt;- c(&quot;SP&quot;, &quot;SC&quot;, &quot;DF&quot;, &quot;MG&quot;, &quot;RS&quot;, &quot;MS&quot;, &quot;PR&quot;, &quot;AC&quot;, &quot;PI&quot;,&quot;TO&quot;, &quot;MT&quot;, &quot;RO&quot;, &quot;GO&quot;, &quot;RJ&quot;, &quot;ES&quot;, &quot;MA&quot;, &quot;PB&quot;, &quot;AL&quot;, &quot;RN&quot;, &quot;PA&quot;, &quot;SE&quot;, &quot;PE&quot;, &quot;RR&quot;, &quot;CE&quot;, &quot;AM&quot;, &quot;BA&quot;, &quot;AP&quot;)
  2. VID &lt;- c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8)
  3. VIDBrazil &lt;- data.frame(State, VID)
  5. spdf_Brasil &lt;- geojson_read(&quot;C:/Users/Giovanni/Desktop/R/brazil_geo.json&quot;, what = &quot;sp&quot;) #Import JSON file
  7. spdf_Brasil_A&lt;- tidy(spdf_Brasil) #Convert to data frame
  9. ggplot() + #plot the map
  10.   geom_polygon(data = spdf_Brasil_A, aes( x = long, y = lat, group = group), fill=&quot;#69b3a2&quot;, color=&quot;white&quot;) +
  11.   theme_void() +
  12.   coord_map()
  14. spdf_Brasil_MVI &lt;- spdf_Brasil_A %&gt;% #add the values of deaths to use for color coding the map
  15.   left_join(. , VIDBrazil, by=c(&quot;id&quot;=&quot;State&quot;))
  18. ggplot() + #ploting with color coding
  19.   geom_polygon(data = spdf_Brasil_MVI, aes( x = long, y = lat, group = group, fill= VID), color=&quot;green&quot;) + 
  20.   theme_void() +
  21.   scale_fill_gradient(name = &quot;Mortes por 100 mil habitantes&quot;, low = &quot;#fdc0b2&quot;, high = &quot;#9a1f03&quot;, limits = c(0, 60.0)) +
  22.   coord_map()

Obtained map

What am I missing? What should I do to get the map color coded?

Thank you for your assistance!

答案1

得分: 0

在broom::tidy()中，GeoJSON中的id字段未被处理，因此在生成的数据框中没有State代码，联接将无法按预期工作。尽管名为id的列仍然存在，因此left_join()不会引发错误。如果检查spdf_Brasil_MVI，您可能会注意到所有VID值都是NA。而id列包含存储为字符的数值，而不是State代码：

  1. glimpse(spdf_Brasil_MVI)
  2. #&gt; 行: 581,207
  3. #&gt; 列: 8
  4. #&gt; $ long  &lt;dbl&gt; -73.33251, -73.27482, -72.87016, -72.69151, -72.65935, -72.18967, -72.17318, -…
  5. #&gt; $ lat   &lt;dbl&gt; -7.324879, -7.350334, -7.528882, -7.610621, -7.624979, -7.720489, -7.723841, -…
  6. #&gt; ...
  7. #&gt; $ id    &lt;chr&gt; &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;…
  8. #&gt; $ VID   &lt;dbl&gt; NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…

使用sf包处理空间数据，并使用ggplot的geom_sf()图层绘制sf对象可以简化这个问题：

  1. library(sf)
  2. library(ggplot2)
  3. library(dplyr)
  4. VIDBrazil <- data.frame(
  5.   State = c("SP", "SC", "DF", "MG", "RS", "MS", "PR", "AC", "PI","TO", "MT", "RO", "GO", "RJ", "ES", "MA", "PB", "AL", "RN", "PA", "SE", "PE", "RR", "CE", "AM", "BA", "AP"),
  6.   VID = c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8))
  8. sf_Brasil <- st_read("brazil_geo.json", as_tibble = TRUE, quiet = TRUE)
  9. # sf对象，包含所有要素（id，name）+ geometry列和投影细节：
  10. sf_Brasil
  11. #&gt; 简单特征集合，27个要素和2个字段
  12. #&gt; 几何类型: MULTIPOLYGON
  13. #&gt; 尺寸:     XY
  14. #&gt; 边界框:  xmin: -73.98971 ymin: -33.74708 xmax: -28.84694 ymax: 5.264878
  15. #&gt; 大地测量参考:  WGS 84
  16. #&gt; # A tibble: 27 &#215; 3
  17. #&gt;    id    name                                                           geometry
  18. #&gt;    &lt;chr&gt; &lt;chr&gt;                                                &lt;MULTIPOLYGON [&#176;]&gt;
  19. #&gt;  1 AC    Acre             (((-73.33251 -7.324879, -73.27482 -7.350334, -72.8701…
  20. #&gt;  2 AL    Alagoas          (((-35.90153 -9.861805, -35.90153 -9.862083, -35.9012…
  21. #&gt;  3 AP    Amap&#225;            (((-50.02403 0.859862, -50.02403 0.859306, -50.02347 …
  22. # ...
  24. p1 <- ggplot(sf_Brasil) +
  25.   geom_sf(fill="#69b3a2", color="white") +
  26.   theme_void()
  28. sf_Brasil_MVI <- left_join(sf_Brasil, VIDBrazil, by = c("id" = "State"))
  29. # join result, includes VID values:
  30. sf_Brasil_MVI
  31. #&gt; 简单特征集合，27个要素和3个字段
  32. #&gt; 几何类型: MULTIPOLYGON
  33. #&gt; 尺寸:     XY
  34. #&gt; 边界框:  xmin: -73.98971 ymin: -33.74708 xmax: -28.84694 ymax: 5.264878
  35. #&gt; 大地测量参考:  WGS 84
  36. #&gt; # A tibble: 27 &#215; 4
  37. #&gt;    id    name                                                     geometry   VID
  38. #&gt;    &lt;chr&gt; &lt;chr&gt;                                          &lt;MULTIPOLYGON [&#176;]&gt; &lt;dbl&gt;
  39. #&gt;  1 AC    Acre             (((-73.33251 -7.324879, -73.27482 -7.350334, -7…  21.2
  40. #&gt;  2 AL    Alagoas          (((-35.90153 -9.861805, -35.90153 -9.862083, -3…  31.8
  41. #&gt;  3 AP    Amap&#225;            (((-50.02403 0.859862, -50.02403 0.859306, -50.…  53.8
  42. # ...
  44. p2 <- ggplot(sf_Brasil_MVI) + 
  45.   geom_sf(aes(fill = VID), color="green") +
  46.   scale_fill_gradient(name = "每十万人口的死亡人数", low = "#fdc0b2", high = "#9a1f03", limits = c(0, 60.0)) +
  47.   theme
  49. <details>
  50. <summary>英文:</summary>
  52. `id` field in GeoJSON is not handled by `broom::tidy()` as-is, so there are no State codes in the resulting data frame and join will not work as intended. Though the column named `id` is still there, thus `left_join()` will not throw an error.   If you check your `spdf_Brasil_MVI`, you may notice that all `VID` values are `NA`s. And instead of State codes, `id` column includes numeric values stored as characters:
  53. ``` r
  54. glimpse(spdf_Brasil_MVI)
  55. #&gt; Rows: 581,207
  56. #&gt; Columns: 8
  57. #&gt; $ long  &lt;dbl&gt; -73.33251, -73.27482, -72.87016, -72.69151, -72.65935, -72.18967, -72.17318, -…
  58. #&gt; $ lat   &lt;dbl&gt; -7.324879, -7.350334, -7.528882, -7.610621, -7.624979, -7.720489, -7.723841, -…
  59. #&gt; ...
  60. #&gt; $ id    &lt;chr&gt; &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;, &quot;1&quot;…
  61. #&gt; $ VID   &lt;dbl&gt; NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…

Using sf package for handling spatial data and ggplot's geom_sf() layer for plotting sf objects simplifies this quite a bit:

  1. library(sf)
  2. library(ggplot2)
  3. library(dplyr)
  4. VIDBrazil &lt;- data.frame(
  5.   State = c(&quot;SP&quot;, &quot;SC&quot;, &quot;DF&quot;, &quot;MG&quot;, &quot;RS&quot;, &quot;MS&quot;, &quot;PR&quot;, &quot;AC&quot;, &quot;PI&quot;,&quot;TO&quot;, &quot;MT&quot;, &quot;RO&quot;, &quot;GO&quot;, &quot;RJ&quot;, &quot;ES&quot;, &quot;MA&quot;, &quot;PB&quot;, &quot;AL&quot;, &quot;RN&quot;, &quot;PA&quot;, &quot;SE&quot;, &quot;PE&quot;, &quot;RR&quot;, &quot;CE&quot;, &quot;AM&quot;, &quot;BA&quot;, &quot;AP&quot;),
  6.   VID = c(7.9, 10.1, 11.2, 11.4, 15.9, 20.7, 20.8, 21.2, 23.8, 24.3, 24.9, 25.0, 26.1, 27.2, 28.2, 28.3, 28.6, 31.8, 32.4, 32.8, 33.9, 34.8, 35.5, 37.0, 39.1, 44.9, 53.8))
  8. sf_Brasil &lt;- st_read(&quot;brazil_geo.json&quot;, as_tibble = TRUE, quiet = TRUE)
  9. # sf object, a data.frame / tibble with all features (id, name) + geometry column
  10. # and projection details:
  11. sf_Brasil
  12. #&gt; Simple feature collection with 27 features and 2 fields
  13. #&gt; Geometry type: MULTIPOLYGON
  14. #&gt; Dimension:     XY
  15. #&gt; Bounding box:  xmin: -73.98971 ymin: -33.74708 xmax: -28.84694 ymax: 5.264878
  16. #&gt; Geodetic CRS:  WGS 84
  17. #&gt; # A tibble: 27 &#215; 3
  18. #&gt;    id    name                                                           geometry
  19. #&gt;    &lt;chr&gt; &lt;chr&gt;                                                &lt;MULTIPOLYGON [&#176;]&gt;
  20. #&gt;  1 AC    Acre             (((-73.33251 -7.324879, -73.27482 -7.350334, -72.8701…
  21. #&gt;  2 AL    Alagoas          (((-35.90153 -9.861805, -35.90153 -9.862083, -35.9012…
  22. #&gt;  3 AP    Amap&#225;            (((-50.02403 0.859862, -50.02403 0.859306, -50.02347 …
  23. # ...
  25. p1 &lt;- ggplot(sf_Brasil) +
  26.   geom_sf(fill=&quot;#69b3a2&quot;, color=&quot;white&quot;) +
  27.   theme_void()
  29. sf_Brasil_MVI &lt;- left_join(sf_Brasil, VIDBrazil, by = join_by(id == State))
  30. # join result, includes VID values:
  31. sf_Brasil_MVI
  32. #&gt; Simple feature collection with 27 features and 3 fields
  33. #&gt; Geometry type: MULTIPOLYGON
  34. #&gt; Dimension:     XY
  35. #&gt; Bounding box:  xmin: -73.98971 ymin: -33.74708 xmax: -28.84694 ymax: 5.264878
  36. #&gt; Geodetic CRS:  WGS 84
  37. #&gt; # A tibble: 27 &#215; 4
  38. #&gt;    id    name                                                     geometry   VID
  39. #&gt;    &lt;chr&gt; &lt;chr&gt;                                          &lt;MULTIPOLYGON [&#176;]&gt; &lt;dbl&gt;
  40. #&gt;  1 AC    Acre             (((-73.33251 -7.324879, -73.27482 -7.350334, -7…  21.2
  41. #&gt;  2 AL    Alagoas          (((-35.90153 -9.861805, -35.90153 -9.862083, -3…  31.8
  42. #&gt;  3 AP    Amap&#225;            (((-50.02403 0.859862, -50.02403 0.859306, -50.…  53.8
  43. # ...
  45. p2 &lt;- ggplot(sf_Brasil_MVI) + 
  46.   geom_sf(aes(fill = VID), color=&quot;green&quot;) +
  47.   scale_fill_gradient(name = &quot;Mortes por 100 mil habitantes&quot;, low = &quot;#fdc0b2&quot;, high = &quot;#9a1f03&quot;, limits = c(0, 60.0)) +
  48.   theme_void() +
  49.   theme(legend.position = &quot;bottom&quot;)
  51. p1;p2

用ggplot制作Cloropleth<!-- -->用ggplot制作Cloropleth<!-- -->

For labelling sf objects we could use geom_sf_label() or geom_sf_text(), but when things can get bit crowded and it's not so easy to fit all labels, geom_label_repel() & geom_text_repel() form ggrepel package are worth considering.

  1. library(ggrepel)
  3. # add layer to existing p2 plot
  4. # get coordinates for geom_label_repel from sf_coordinates stats
  5. # name : column in p2 data (sf_Brasil_MVI) 
  6. # geometry : geometry column of p2 data (sf_Brasil_MVI)
  7. p2 + geom_label_repel(aes(label = name, geometry = geometry),
  8.                       stat = &quot;sf_coordinates&quot;)

用ggplot制作Cloropleth<!-- -->

<sup>Created on 2023-03-12 with reprex v2.0.2</sup>

huangapple
  • 本文由 发表于 2023年3月12日 14:25:05
  • 转载请务必保留本文链接：https://go.coder-hub.com/75711401.html
  • ggplot2
  • maps
  • r
  • visualization
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定