How to overlay a 2d density plot on top of a map-CodePudding

I need some help in overlaying a 2d density plot on top of a ggmap plot. I don't really know how to procede. Any help is welcome.

This is a subset of my data, including start and end coordinates for each individual ride of a fictional bike sharing company:

df <- structure(list(start_lat = c(41.94018, 41.890762, 41.845695, 
41.857813, 41.9287386666667, 42.0044803333333, 41.879255, 41.886835, 
41.874734, 41.95469, 41.95, 41.8809518333333, 41.96590013976, 
41.909668, 41.931248, 41.96167, 41.912133, 41.87947235235, 41.936266, 
41.922695, 41.9101756666667, 41.86, 41.91468, 41.892278, 42.03, 
41.911386, 41.9716, 41.93, 41.940195, 41.9560855, 41.915784, 
41.93314, 41.943739, 41.8671848333333, 41.87464, 41.882242, 41.926277, 
41.96167, 41.76, 41.883668, 41.967096, 41.8, 41.9024035, 41.939743, 
41.9093960065, 41.915983, 41.87772613, 41.8984238333333, 41.8836331666667, 
41.925905, 41.967096, 41.92, 41.884576228, 41.838499, 41.9028846666667, 
41.89993001, 41.8, 41.866095, 41.97, 41.9093960065, 41.88, 41.8922376666667, 
41.81, 41.9438251666667, 41.883668, 41.9207793333333, 41.954383, 
41.9434726666667, 41.8945555, 41.911386, 41.88917683258, 41.86722595682, 
41.8531223333333, 41.92, 41.919936, 41.90096, 41.894722, 41.872187, 
41.881892, 41.920082, 41.897448, 41.88917683258, 41.9, 41.925858, 
41.89, 41.8908470406238, 41.85, 41.890173, 41.92556258, 41.885637, 
41.9030376666667, 41.93314, 41.838198, 41.892278, 41.93, 41.894722, 
41.90345, 41.6922943333333, 41.9080621666667, 42.025784), start_lng = c(-87.65304, 
-87.631697, -87.6225141666667, -87.62455, -87.6538015, -87.6615086666667, 
-87.639904, -87.62232, -87.6498425, -87.67393, -87.71, -87.6167566666667, 
-87.6936384935, -87.648128, -87.644336, -87.65464, -87.634656, 
-87.6256886059, -87.652662, -87.697153, -87.6823075, -87.62, 
-87.64332, -87.612043, -87.71, -87.638677, -87.650154, -87.74, 
-87.6529666666667, -87.668857, -87.634581, -87.64776, -87.66402, 
-87.6260033333333, -87.65703, -87.641066, -87.630834, -87.65464, 
-87.58, -87.64867, -87.667429, -87.58, -87.6277486666667, -87.658865, 
-87.6776919292, -87.677335, -87.65478743, -87.6223878333333, 
-87.629143, -87.64926, -87.667429, -87.7, -87.63188991, -87.6080766666667, 
-87.6874035, -87.63443007, -87.59, -87.607267, -87.71, -87.6776919292, 
-87.63, -87.6119485, -87.61, -87.671138, -87.64867, -87.6637163333333, 
-87.648043, -87.6796343333333, -87.6534645, -87.638677, -87.6385057718, 
-87.6153553902, -87.6318963333333, -87.74, -87.64883, -87.623777, 
-87.634362, -87.661501, -87.648789, -87.677855, -87.628722, -87.6385057718, 
-87.62, -87.638973, -87.66, -87.6186168193817, -87.72, -87.626185, 
-87.65840426, -87.641823, -87.631299, -87.64776, -87.645143, 
-87.612043, -87.71, -87.634362, -87.667747, -87.6426485, -87.6315093333333, 
-87.684107), end_lat = c(41.918306, 41.886875, 41.8456825, 41.8530845574128, 
41.890831, 41.99, 41.885637, 41.881319815, 41.88, 41.961068, 
41.93, 41.880958, 41.966399801841, 41.89766, 41.9267559875, 41.9578665241517, 
41.911386, 41.867888, 41.95078, 41.932588, 41.9245285, 41.8776751666667, 
41.9105780349, 41.9239313113662, 42.0192226666667, 41.904613, 
41.9947796884, 41.96, 41.9296915, 41.94, 41.94, 41.907066, 41.923931, 
41.8707831666667, 41.87772613, 41.872187, 41.892278, 41.961004, 
41.7689161666667, 41.8793563587, 41.95078, 41.79, 41.882242, 
41.932225, 41.912133, 41.9093960065, 41.8810317, 41.9, 41.89, 
41.912133, 41.926277, 41.93190196886, 41.874053, 41.8368228333333, 
41.882754, 41.894666, 41.76, 41.882134, 41.96, 41.89637337, 41.87, 
41.8787191666667, 41.8, 41.91, 41.917805, 41.88, 41.926277, 41.93, 
41.8990156666667, 41.890762, 41.8854833079, 41.874754, 41.85, 
41.94, 41.920771, 41.894345, 41.94334, 41.871737, 41.88338, 41.92154, 
41.882134, 41.902997, 41.876243, 41.892278, 41.89, 41.886024, 
41.86, 41.8918473721099, 41.9093960065, 41.8854833079, 41.89, 
41.92883, 41.834734, 41.891466, 41.9296816666667, 41.902973, 
41.918491153687, 41.75, 41.9218326666667, 41.9840446107), end_lng = c(-87.636282, 
-87.62603, -87.6224476666667, -87.6319313049316, -87.6313945, 
-87.66, -87.641823, -87.6295209193, -87.65, -87.695439, -87.71, 
-87.616743, -87.6887042820454, -87.62351, -87.6344287848, -87.6495051383972, 
-87.638677, -87.623041, -87.659172, -87.636427, -87.658447, -87.6240391666667, 
-87.6494219288, -87.6358245313168, -87.6736431666667, -87.640552, 
-87.6602845349, -87.69, -87.7080808333333, -87.67, -87.68, -87.667252, 
-87.635825, -87.6257745, -87.65478743, -87.661501, -87.612043, 
-87.649603, -87.634775, -87.6297910363, -87.659172, -87.6, -87.641066, 
-87.658617, -87.634656, -87.6776919292, -87.62408432, -87.62, 
-87.63, -87.634656, -87.630834, -87.7011951301, -87.627716, -87.6133453333333, 
-87.6259215, -87.638437, -87.55, -87.625125, -87.69, -87.66098386, 
-87.62, -87.6355345, -87.59, -87.66, -87.682437, -87.63, -87.630834, 
-87.71, -87.6299358333333, -87.631697, -87.6523048564, -87.649807, 
-87.64, -87.73, -87.663712, -87.622798, -87.67097, -87.65103, 
-87.64117, -87.653818, -87.625125, -87.683825, -87.624426, -87.612043, 
-87.65, -87.624117, -87.72, -87.6205801963806, -87.6776919292, 
-87.6523048564, -87.63, -87.668507, -87.625813, -87.626761, -87.7081071666667, 
-87.63128, -87.6974228024483, -87.64, -87.6439593333333, -87.6602738295
)), row.names = c(NA, -100L), class = "data.frame")

Next, i extracted min and max values for latitude and longitude and then used those values as limits for my map plot:

library(ggmap)
library(ggplot2)

map_lim <- data.frame(
  min_lat = min(df[, c("start_lat", "end_lat")]), 
  max_lat = max(df[, c("start_lat", "end_lat")]),
  min_lng = min(df[, c("start_lng", "end_lng")]), 
  max_lng = max(df[, c("start_lng", "end_lng")])
)

map_lim

##    min_lat max_lat min_lng max_lng
## 1 41.69229   42.03  -87.74  -87.55

map <- get_stamenmap(
  bbox = c(left = map_lim$min_lng, right = map_lim$max_lng, bottom = map_lim$min_lat, top = map_lim$max_lat)
)

ggmap(map)

This is the resulting map:

Now i want to add a 2d density plot based on the same set of coordinates on top of this map, but i don't know how to do it. This is my 2d density plot code.

density2d <- ggplot(df, aes(x = start_lng, y = start_lat))  
  coord_equal(xlim = c(map_lim$min_lng, map_lim$max_lng), ylim = c(map_lim$min_lat, map_lim$max_lat))  
  xlab("Longitude")  
  ylab("Latitude")  
  stat_density2d(aes(fill = ..level..), alpha = 0.5, geom = "polygon")  
  scale_fill_viridis_c()
  
density2d

CodePudding user response：

You simply add the density layer to your map:

ggmap(map)  
  stat_density2d(data = df, aes(x = start_lng, y = start_lat,
                 fill = ..level..), alpha = 0.5, geom = "polygon")  
  scale_fill_viridis_c()