Подсчет плотности в тепловых картах

ИСПРАВЛЕНО

Я понял, что мой предыдущий ответ нужно пересмотреть. Итак, вот оно. Если вы хотите узнать, сколько точек данных существует на каждом уровне контура, на самом деле вам нужно сделать много вещей. Если вы согласны использовать опцию leaflet ниже, ваша жизнь станет намного проще.

Во-первых, давайте получим карту Детройта и создадим образец фрейма данных.

library(dplyr)
library(ggplot2)
library(ggmap)

mymap <- get_map(location = "Detroit", zoom = 8)

### Create a sample data
set.seed(123)
mydata <- data.frame(long = runif(min = -84, max = -82.5, n = 100),
                     lat = runif(min = 42, max = 42.7, n = 100))

Теперь мы рисуем карту и сохраняем ее как g.

g <- ggmap(mymap) +
     stat_density2d(data = mydata,
                    aes(x = long, y = lat, fill = ..level..),
                    size = 0.5, bins = 10, geom = "polygon")

Настоящая работа начинается здесь. Чтобы узнать количество точек данных на всех уровнях, вы хотите использовать фрейм данных, который генерирует ggplot. В этом фрейме данных у вас есть данные для полигонов. Эти многоугольники используются для рисования линий уровня. Вы можете видеть это на следующем изображении, где я рисую три уровня на карте.

### Create a data frame so that we can find how many data points exist
### in each level.

mydf <- ggplot_build(g)$data[[4]]

### Check where the polygon lines are. This is just for a check.

check <- ggmap(mymap) +
         geom_point(data = mydata, aes(x = long, y = lat)) +
         geom_path(data = subset(mydf, group == "1-008"), aes(x = x, y = y)) +
         geom_path(data = subset(mydf, group == "1-009"), aes(x = x, y = y)) +
         geom_path(data = subset(mydf, group == "1-010"), aes(x = x, y = y)) 

Следующим шагом является создание вектора уровня для легенды. Мы группируем данные по группам (например, 1-010) и берем первую строку для каждой группы, используя slice(). Затем разгруппируйте данные и выберите второй столбец. Наконец, создайте вектор с unlist(). Мы возвращаемся к lev в конце.

mydf %>%
group_by(group) %>%
slice(1) %>%
ungroup %>%
select(2) %>%
unlist -> lev

Теперь мы разделяем данные полигона (т. е. mydf) по группам и создаем полигон для каждого уровня. Поскольку у нас 11 уровней (11 полигонов), мы используем lapply(). В нахлестной петле нам нужно сделать; 1) извлечь столбец для долготы и широты, 2) создать полигон, 3) преобразовать полигоны в пространственные полигоны, 4) назначить CRS, 5) создать фиктивный фрейм данных и 6) создать SpatialPolygonsDataFrames.

mylist <- split(mydf, f = mydf$group)

test <- lapply(mylist, function(x){

              xy <- x[, c(3,4)]

              circle <- Polygon(xy, hole = as.logical(NA))

              SP <- SpatialPolygons(list(Polygons(list(circle), ID = "1")))

              proj4string(SP) <- CRS("+proj=longlat +ellps=WGS84")

              df <- data.frame(value = 1, row.names = "1")

              circleDF <- SpatialPolygonsDataFrame(SP, data = df)

            })

Теперь вернемся к исходным данным. Нам нужно преобразовать фрейм данных в SpatialPointsDataFrame. Это потому, что нам нужно разделить данные и найти, сколько точек данных существует в каждом полигоне (на каждом уровне). Во-первых, получите долготу и широту из вашего data.frame. Убедитесь, что заказ указан в долготе/широте.

xy <- mydata[,c(1,2)]

Затем мы создаем SPDF (SpatialPolygonsDataFrame). Вы хотите иметь идентичную строку proj4string между пространственными полигонами и данными пространственных точек.

spdf <- SpatialPointsDataFrame(coords = xy, data = mydata,
                               proj4string = CRS("+proj=longlat +ellps=WGS84"))

Затем мы подмножаем данные (mydata), используя каждый полигон.

ana <- lapply(test, function(y){

              mydf <- as.data.frame(spdf[y, ])

            })

Точки данных перекрываются на разных уровнях; у нас дублирование. Сначала мы пытаемся найти уникальные точки данных для каждого уровня. Мы связываем фреймы данных в ana и создаем фрейм данных foo1. Мы также создаем фрейм данных, в котором мы хотим найти уникальное количество точек данных. Мы следим за тем, чтобы имена столбцов были одинаковыми между foo1 и foo2. Используя setdiff() и nrow(), мы можем найти уникальное количество точек данных на каждом уровне.

total <- lapply(11:2, function(x){

                foo1 <- bind_rows(ana[c(11:x)])
                foo2 <- as.data.frame(ana[x-1])
                names(foo2) <- names(foo1)
                nrow(setdiff(foo2, foo1))               
              })

Наконец, нам нужно найти количество точек данных для самого внутреннего уровня, то есть уровня 11. Мы выбираем фрейм данных для уровня 11 в ana, создаем фрейм данных и подсчитываем количество строк.

 bob <- nrow(as.data.frame(ana[11]))
 out <- c(bob,unlist(total))

 ### check if total is 100
 ### sum(out)
 ### [1] 100

Мы назначаем перевернутые out в качестве имен для lev. Это потому, что мы хотим показать, сколько точек данных существует для каждого уровня в легенде.

 names(lev) <- rev(out)

Теперь мы готовы добавить легенду.

 final <- g +
          scale_fill_continuous(name = "Total",
                                guide = guide_legend(),
                                breaks = lev)

 final

ВАРИАНТ ЛИСТКИ

Если вы используете пакет листовок, вы можете группировать точки данных с разным увеличением. Листовка подсчитывает точки данных в определенных областях и указывает числа в кружках, как показано на следующем рисунке. Чем больше вы увеличиваете масштаб, тем больше листовка разбивает точки данных на небольшие группы. С точки зрения рабочей нагрузки, это намного легче. Кроме того, ваша карта интерактивна. Это может быть лучшим вариантом.

library(leaflet)
leaflet(mydf) %>%
addTiles() %>%
addMarkers(clusterOptions = markerClusterOptions())