如何绘制每个多边形内的观测结果图？

Question

我想用下面的多边形绘制一张choropleth图：

library(sp)
library(sf)
library(spatialEco)
library(tigris)

br_tracts <- tracts(state = 'LA', county = 'East Baton Rouge', cb = T, year = 2016)
plot(br_tracts)

​

​

然后，我想将颜色映射到每个多边形内的观察点的数量。以下是我的示例数据。请注意，您需要上传tigris包才能创建dtSpatial。

dtSpatial <- new("SpatialPointsDataFrame"
    , data = structure(list(parenttype = c("garbage", "garbage", "garbage", 
"garbage", "garbage", "garbage", "garbage", "garbage", "garbage", 
"garbage", "garbage", "garbage", "garbage", "garbage", "garbage", 
"garbage", "garbage", "garbage", "garbage", "garbage")), .Names = "parenttype", row.names = c(NA, 
-20L), class = c("tbl_df", "tbl", "data.frame"))
    , coords.nrs = numeric(0)
    , coords = structure(c(-91.043777, -91.026382, 0, -91.027748, 0, -91.08049, 
-91.047173, -91.172501, -91.162384, -91.139465, -91.087585, -91.152748, 
-91.163086, -91.185814, -91.135101, 0, -91.105972, 0, -91.168846, 
-91.041435, 30.452148, 30.447191, 0, 30.412008, 0, 30.415289, 
30.420155, 30.430065, 30.478041, 30.460482, 30.429127, 30.469275, 
30.436682, 30.420218, 30.453447, 0, 30.431898, 0, 30.466148, 
30.416723), .Dim = c(20L, 2L), .Dimnames = list(NULL, c("longitude", 
"latitude")))
    , bbox = structure(c(-91.185814, 0, 0, 30.478041), .Dim = c(2L, 2L), .Dimnames = list(
    c("longitude", "latitude"), c("min", "max")))
    , proj4string = new("CRS"
    , projargs = NA_character_
)
)

为了查看给定的观察值是否出现在多边形中，我尝试：

pts.poly <- point.in.polygon(point.x = dtSpatial$longitude, 
                          point.y = dtSpatial$latitude,
                          pol.x = fortify(br_tracts)$long,
                          pol.y = fortify(br_tracts)$lat )

> pts.poly
 [1] 1 1 0 0 0 0 0 1 1 0 1 0 0 0 0 0 1 0 0 0

但是如何将颜色映射到每个多边形中的观测值数量？

Answer 1

这里有一种适合你的方法。看到数据后，您要做的第一件事就是向dtSpatial添加一个投影。然后，您需要计算每个多边形中有多少个数据点。您可以使用GISTools包中的poly.counts()实现这一点。输出是一个向量。由于该函数使用多边形ID来计算数据点，因此您将看到ID作为名称。您希望使用stack()将其转换为数据框。您可以将br_tracts转换为ggplot的数据框。最后，绘制一张地图。我使用了geom_cartogram()。我使用了这个函数两次。我先画多边形。然后，我使用count用颜色填充它们。我正在匹配map_id来完成这项工作。注意:您可以使用geom_map()、geom_polygon()。coord_map()和theme_map()是可选的。

library(tigris)
library(GISTools)
library(RColorBrewer)
library(ggalt)
library(ggthemes)

# Add projection to dtSpatial, which is identical to projection of br_tracts
proj4string(dtSpatial) <- CRS(proj4string(br_tracts))

# How many data poiints stay in each polygon. Polygon ID is the name of the vector.

count <- poly.counts(pts = dtSpatial, polys = br_tracts)
count <- stack(count)

mymap <- fortify(br_tracts)

ggplot() +
geom_cartogram(data = mymap, aes(x = long, y = lat, map_id = id), 
               map = mymap) +
geom_cartogram(data = count, aes(fill = values, map_id = ind),
               map = mymap, color = "black", size = 0.3) +
scale_fill_gradientn(colours = rev(brewer.pal(10, "Spectral"))) +
coord_map() +
theme_map()

​

​