Search code examples
rggplot2gisggmapdensity-plot

Density Heatmap with 1 Dimension Density Plot Alignment


I am using the method outlined in this post in order to create a density heatmap on a map paired with the one dimensional x and y density plots. I was able to produce the following image: Density Heatmap

However, I am having difficulty aligning the shaded regions in the top density plot with the heatmap on the map. I discovered that this has to to with the fact that I am applying a mercator projection in the code to produce the heatmap, whereas the top and right density plots are just using the coordinates from the data.

Is there a way that I can align the projected coordinates in the heatmap with the nonprojected coordinates in the one dimensional density plot?

Below is the relevant section of my code (I used the grid/gtable packages to combine the different graphs).

EDIT: The data is comprised of three columns: UniqueID, Lon, and Lat. In this case all the Lon/Lat pairs are in/around New York City, but for troubleshooting purposes the exact pairs don't matter as long as they can be projected.

library(gtable)
library(ggplot2)   
library(ggmap) 

ny <- get_map(location = "Washington Heights, New York City", zoom = 12, color = "bw")

    ggm <- ggmap(ny, extent = "normal", maprange = F) %+% data +
          aes(x = Lon, y = Lat) +
          stat_density2d(aes(fill = ..level..,
                             alpha = ..level..), color = "darkgreen",
                         geom = "polygon", show.legend = TRUE) +
          geom_polygon(aes(x, y), 
                       data.frame(x = c(Inf, Inf, -73.90, -73.90), y = c(Inf, 40.84, 40.84, Inf)),
                       alpha = 0.5, colour = NA, fill = "red") +
          coord_map(projection = "mercator",
                    xlim = c(attr(ny, "bb")$ll.lon, attr(ny, "bb")$ur.lon),
                    ylim = c(attr(ny, "bb")$ll.lat, attr(ny, "bb")$ur.lat))

    xd <- data.frame(density(data$Lon)[c("x", "y")])
    gg1 <- ggplot(xd, aes(x, y)) +
      theme(panel.grid.minor = element_line(colour = NA),
                panel.background = element_rect(fill = NA, colour = NA)) +
      labs(y = "Density") +
      geom_area(data = subset(xd, x > -73.90), fill = "red") +
      geom_line() +
      coord_cartesian(c(attr(ny, "bb")$ll.lon, attr(ny, "bb")$ur.lon))

image <- gtable_filter(ggplotGrob(ggm), pattern = "panel", trim = TRUE, fixed=TRUE)
image <- gtable_add_cols(image, unit(0.2, "null"), 1)
image <- gtable_add_grob(image, gtable_filter(ggplotGrob(gg1), pattern = "panel", trim = TRUE, fixed=TRUE), 1, 1) 

Solution

  • I don't think the problem is to do with the projection. If you check the x-range for the top density plot and the x-range for the map, I don't think they are the same. Similarly, the y-range for the right density plot compared to the y-range of the map. See to the bottom for code. The reason the ranges are not the same is that, by default, ggmap has no expansion factor for the axes whereas ggplot, by default, has an expansion factor. Set expand = c(0,0) in the two density plots.

    There seems to be another issue. The right density plot needs to be flipped, but the flipping interferes with the limits set via coord_cartesian. If I set the limits using xlim (and ylim), then the limits remain after the flipping.

    library(ggmap)
    library(gtable)
    library(ggplot2)    
    library(grid)
    
    
    ny <- get_map(location = "Washington Heights, New York City", zoom = 12, color = "bw")
    
    # Pretend data
    set.seed(5126)
    df = data.frame(Lon = rnorm(100, mean = -73.90, sd = .015), Lat = rnorm(100, mean = 40.84, sd = .015))
    
    # Map
       ggm <- ggmap(ny, extent = "normal", maprange = F) %+% df +
              aes(x = Lon, y = Lat) +  
              stat_density2d(aes(fill = ..level..,
                                 alpha = ..level..),
                             geom = "polygon", show.legend = TRUE) +
              geom_polygon(aes(x, y), 
                           data.frame(x = c(Inf, Inf, -73.90, -73.90), y = c(Inf, 40.84, 40.84, Inf)),
                           alpha = 0.5, colour = NA, fill = "red") + geom_point() +
              coord_map(projection = "mercator", 
                        xlim = c(attr(ny, "bb")$ll.lon, attr(ny, "bb")$ur.lon),
                        ylim = c(attr(ny, "bb")$ll.lat, attr(ny, "bb")$ur.lat)) 
    
    
    # Top density plot
    xd <- data.frame(density(df$Lon)[c("x", "y")])
    
    gg1 <- ggplot(xd, aes(x, y)) +
      theme(panel.grid.minor = element_line(colour = NA),
            panel.background = element_rect(fill = NA, colour = NA)) +
      labs(y = "Density") +
      geom_area(data = subset(xd, x > -73.90), fill = "red") +
      geom_line() +
      scale_x_continuous(limits = c(attr(ny, "bb")$ll.lon, attr(ny, "bb")$ur.lon), expand = c(0,0)) 
    
    image <- gtable_filter(ggplotGrob(ggm), pattern = "panel", trim = TRUE, fixed=TRUE)
    image <- gtable_add_rows(image, unit(0.2, "null"), 0)
    image <- gtable_add_grob(image, gtable_filter(ggplotGrob(gg1), pattern = "panel", trim = TRUE, fixed=TRUE), 1, 1) 
    
    grid.newpage()
    grid.draw(image)
    
    
    # Right density plot
    xd <- data.frame(density(df$Lat)[c("x", "y")])
    
    gg2 <- ggplot(xd, aes(x, y)) +
       theme(panel.grid.minor = element_line(colour = NA),
             panel.background = element_rect(fill = NA, colour = NA)) +
       labs(y = "Density") +
       geom_area(data = subset(xd, x > 40.84), fill = "red") + 
       geom_line() + 
       scale_x_continuous(limits = c(attr(ny, "bb")$ll.lat, attr(ny, "bb")$ur.lat), expand = c(0,0)) +
       coord_flip() 
    
    image <- gtable_add_cols(image, unit(0.2, "null"), 1)
    image <- gtable_add_grob(image, gtable_filter(ggplotGrob(gg2), pattern = "panel", trim = TRUE, fixed=TRUE), 2, 2) 
    
    grid.newpage()
    grid.draw(image)
    

    Edit Updating to ggplot2 ver 3.0.0

    # Check that the ranges for the density plots are the same as the ranges for the map
    top = ggplot_build(gg1)
    right = ggplot_build(gg2)
    map = ggplot_build(ggm)
    
    top$layout$panel_params[[1]]$x.range
    map$layout$panel_params[[1]]$x.range
    
    right$layout$panel_params[[1]]$y.range
    map$layout$panel_params[[1]]$y.range
    

    enter image description here