Wie choropleth von Beobachtungen innerhalb jedes Polygons?

Question

Ich möchte eine Choroplethenkarte plotten folgende Polygone:

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

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

Dann würde ich Farbe auf die Anzahl der beobachteten Punkte abbilden möchten, die innerhalb jedes Polygons fallen. Das Folgende sind meine Beispieldaten. Beachten Sie, dass Sie das Tigris-Paket hochladen müssen, um dtSpatial zu erstellen.

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

Um zu sehen, ob eine gegebene Beobachtung in einem Polygon erscheint, habe ich versucht:

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 

Aber wie ordne ich die Farbe auf die Anzahl der Beobachtungen innerhalb jedes Polygons?

Answer 1

Hier ist ein Ansatz für Sie. Wenn Sie Ihre Daten sehen, müssen Sie zuerst eine Projektion zu dtSpatial hinzufügen. Dann möchten Sie zählen, wie viele Datenpunkte in jedem Polygon vorhanden sind. Sie erreichen dies mit poly.counts() im GISTools-Paket. Die Ausgabe ist ein Vektor. Da die Funktion die Polygon-ID zum Zählen von Datenpunkten verwendet, sehen Sie die IDs als Namen. Sie möchten es unter Verwendung von stack() in einen Datenrahmen konvertieren. Sie konvertieren br_tracts in einen Datenrahmen für ggplot. Zum Schluss zeichnen Sie eine Karte. Ich habe geom_cartogram() verwendet. Ich verwende die Funktion zweimal. Ich zeichne zuerst die Polygone. Dann fülle ich sie mit Farben unter Verwendung count. Ich passe map_id an, um diesen Job zu erledigen. Beachten Sie, dass Sie geom_map(), geom_polygon() verwenden können. coord_map() und theme_map() sind 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()