You can plot the map as a scatter3d plot. Create xx and yy arrays with the same dimensions as gg, and with the correct geocoordinates for the map you have retrieved. Then use add_trace(p, x=xx, y=yy, z=0, type=“scatter3d”, color=gg). If you are mapping a large area, gg will be large and you may need to subsample a lower-resolution version to get it to plot correctly.
Here’s my code (I’ve used “landmap” instead of “gg” and “MapLocation” instead of “loc”, and have used “ind1” and “ind2” to reduce the resolution of my map):
landmap <- get_map(location=MapLocation, crop=TRUE)
bb <- as.numeric(attributes(landmap)$bb)
landy <- seq(bb[3], bb[1], length.out = dim(landmap)[1])
landx <- seq(bb[2], bb[4], length.out = dim(landmap)[2])
landxx <- (array(landx, dim=dim(t(landmap))))
landyy <- t(array(landy, dim=dim(landmap)))
ind1 <- seq(1, dim(landxx)[1], 4)
ind2 <- seq(1, dim(landxx)[2], 4)
p <- add_trace(p, x=(c(landxx[ind1, ind2])),
y=(c(landyy[ind1, ind2])),
z=0,
type="scatter3d",
mode="markers",
marker=list(color=c(landmap[ind1, ind2]), size=4),
showlegend=FALSE)