problem using `facet_wrap` in `ggplot` and `ggmap` to show multiple graph

I have a dataset that includes the information shown below:

  Postal.Code   top LATITUDE LONGITUDE
   <chr>       <int>    <dbl>     <dbl>
 1 V6Z 3E6         1     49.27369     -123.130006

structure(list(Postal.Code = c("V6Z 3E6", "V6H 0A7", "V6G 1X5", 
"V6J 2A4", "V6B 2L3", "V5L 4N3", "V6Z 2Y2", "V5V 3N2", "V6P 6Z1", 
"V6E 4R7", "V6B 1S3", "V6Z 3G2", "V6J 0B4", "V5K 0C5", "V5N 5P9", 
"V5L 3S4", "V6Z 2E9", "V6B 1R1", "V6E 1G2", "V6A 1B5", "V6E 4L8", 
"V6E 4L8", "V6E 4L8", "V6B 6A7", "V6B 4K8", "V6E 0B1", "V6Z 3E1", 
"V6E 0B1", "V6Z 1C2", "V6E 0B1", "V6Z 3A4", "V6E 1H1", "V6A 4K7", 
"V6Z 2Y4", "V5T 3G2", "V6E 4R7", "V6E 0B1", "V6Z 3B8", "V6B 1B7", 
"V6B 6M2"), top = c(1L, 4L, 5L, 5L, 3L, 5L, 1L, 5L, 4L, 3L, 4L, 
1L, 3L, 5L, 5L, 3L, 3L, 1L, 5L, 4L, 2L, 3L, 4L, 1L, 1L, 4L, 1L, 
5L, 2L, 4L, 1L, 5L, 1L, 1L, 4L, 2L, 4L, 1L, 5L, 5L), LATITUDE = c(49.27369, 
49.2632107, 49.282456, 49.261313, 49.281036, 49.283819, 49.273184, 
49.254803, 49.2030784, 49.2785751, 49.276925, 49.273926, 49.263131, 
48.004134, 49.284064, 49.272736, 49.279232, 49.280055, 49.2868418, 
49.283923, 49.2786587, 49.2786587, 49.2786587, 49.274144, 49.2750068, 
49.281987, 49.278667, 49.281987, 49.276812, 49.281987, 49.273718, 
49.283478, 49.276, 49.273557, 49.25984, 49.2785751, 49.281987, 
49.273347, 49.2754447, 49.2789398), LONGITUDE = c(-123.130006, 
-123.128055, -123.138561, -123.068251, -123.10882, -123.059851, 
-123.119397, -123.10099, -123.1348227, -123.1302609, -123.118491, 
-123.127924, -123.150686, -113.331951, -123.070319, -123.073661, 
-123.1269454, -123.10569, -123.1301548, -123.101285, -123.1303935, 
-123.1303935, -123.1303935, -123.124925, -123.1257866, -123.125272, 
-123.130037, -123.125272, -123.131891, -123.125272, -123.117243, 
-123.127271, -123.101171, -123.119595, -123.100871, -123.1302609, 
-123.125272, -123.127828, -123.1249486, -123.1158462)), row.names = c(NA, 
-40L), class = c("tbl_df", "tbl", "data.frame"))

the column top is basically a topic number column and I want to plot points on the map for each of topic number. Instead of subsetting the data for each topic, I decided to use facet_wrap as shown below:

van <- get_map(location = "vancouver", zoom = 12)
HoustonMap <- ggmap(van, base_layer = ggplot(aes(x = LONGITUDE, y = LATITUDE),
                                             data = data))
HoustonMap +
  stat_density2d(aes(x = LONGITUDE, y = LATITUDE, fill = stat(level), alpha = ..level..),
                 bins = 10, geom = "polygon",
                 data = data) +
  scale_fill_gradient(low = "black", high = "red") + facet_wrap(~ top)

However, this doesn't show the result properly. the plot is shown below:

It just shows a red point for topic number 2 and that's it. How can I fix this?

Solution

Your code is fine, albeit you have one point beyond the bounds of this map.

The problem is rather the low density of points for each facet that gets mapped to alpha = ..level.. in stat_density2d. To improve visualization, simply raise the number of bins.

library(tidyverse)
library(ggmap)
df_points <- 
  tibble(
    Postal.Code = c("V6Z 3E6", "V6H 0A7", "V6G 1X5", "V6J 2A4", "V6B 2L3", "V5L 4N3", "V6Z 2Y2", "V5V 3N2", "V6P 6Z1", "V6E 4R7", "V6B 1S3", "V6Z 3G2", "V6J 0B4", "V5K 0C5", "V5N 5P9", "V5L 3S4", "V6Z 2E9", "V6B 1R1", "V6E 1G2", "V6A 1B5", "V6E 4L8", "V6E 4L8", "V6E 4L8", "V6B 6A7", "V6B 4K8", "V6E 0B1", "V6Z 3E1", "V6E 0B1", "V6Z 1C2", "V6E 0B1", "V6Z 3A4", "V6E 1H1", "V6A 4K7", "V6Z 2Y4", "V5T 3G2", "V6E 4R7", "V6E 0B1", "V6Z 3B8", "V6B 1B7", "V6B 6M2"), 
    top = c(1L, 4L, 5L, 5L, 3L, 5L, 1L, 5L, 4L, 3L, 4L, 1L, 3L, 5L, 5L, 3L, 3L, 1L, 5L, 4L, 2L, 3L, 4L, 1L, 1L, 4L, 1L, 5L, 2L, 4L, 1L, 5L, 1L, 1L, 4L, 2L, 4L, 1L, 5L, 5L), 
    LATITUDE = c(49.27369, 49.2632107, 49.282456, 49.261313, 49.281036, 49.283819, 49.273184, 49.254803, 49.2030784, 49.2785751, 49.276925, 49.273926, 49.263131, 48.004134, 49.284064, 49.272736, 49.279232, 49.280055, 49.2868418, 49.283923, 49.2786587, 49.2786587, 49.2786587, 49.274144, 49.2750068, 49.281987, 49.278667, 49.281987, 49.276812, 49.281987, 49.273718, 49.283478, 49.276, 49.273557, 49.25984, 49.2785751, 49.281987, 49.273347, 49.2754447, 49.2789398), 
    LONGITUDE = c(-123.130006, -123.128055, -123.138561, -123.068251, -123.10882, -123.059851, -123.119397, -123.10099, -123.1348227, -123.1302609, -123.118491, -123.127924, -123.150686, -113.331951, -123.070319, -123.073661, -123.1269454, -123.10569, -123.1301548, -123.101285, -123.1303935, -123.1303935, -123.1303935, -123.124925, -123.1257866, -123.125272, -123.130037, -123.125272, -123.131891, -123.125272, -123.117243, -123.127271, -123.101171, -123.119595, -123.100871, -123.1302609, -123.125272, -123.127828, -123.1249486, -123.1158462))

vancouver <- get_map(location = "vancouver", zoom = 12)

ggmap(ggmap = vancouver,
      maprange = FALSE, 
      base_layer = ggplot(data = df_points,
                          aes(x = LONGITUDE, y = LATITUDE))) +
  stat_density2d(aes(x = LONGITUDE, y = LATITUDE, fill = ..level.., alpha = ..level..),
                 bins = 1000, # Modify this number to tweak the density areas.
                 geom = "polygon") +
  scale_fill_gradient(low = "red", high = "black") +
  facet_wrap(~top)