How to plot choropleth of observations within each polygon?

Question

I would like to plot a choropleth map using the following polygons:

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

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

Then, I would like to map color to the number of observed points that fall within each polygon. The following is my sample data. Note that you need to upload the tigris package in order to create 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_
)
)

To see whether a given observation appears in a polygon, I tried:

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

But how do I map the color to the number of observations within each polygon?

Answer 1

Here is one approach for you. Seeing your data, the first thing you want to do is to add a projection to dtSpatial. Then, you want to count how many data points exist in each polygon. You can achieve this using poly.counts() in the GISTools package. The output is a vector. Since the function is using polygon ID to count data points, you see the IDs as names. You want to convert it to a data frame using stack(). You convert br_tracts to a data frame for ggplot. Finally, you draw a map. I used geom_cartogram(). I am using the function twice. I am drawing the polygons first. Then, I am filling in them with colors using count. I am matching map_id to do this job. Note you can use geom_map(), geom_polygon(). coord_map() and theme_map() are optional.

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()