Search code examples
rggplot2ggmapgeom-raster

geom_raster() produces a whitish surface ontop of the map


I am trying to plot a heatmap ontop of a geographical map to show the geographic distribution of a variable. The minimum working code, with absurd data, is the following:

library(ggmap)
library(osmdata)
box <- c(left = 2.075, bottom = 41.325, right = 2.25, top = 41.47)
map <- get_stamenmap(bbox = box, maptype = "terrain-lines", zoom = 13)

lon_grid <- seq(2.075, 2.25, length.out = 30)
lat_grid <- seq(41.325, 41.47, length.out = 30)
grid <- expand.grid(lon_grid, lat_grid)
z <- c(rep(NA, 30^2/2), rnorm(30^2/2))
dataset <- cbind(grid, z)

ggmap(map) ### Plot 1

ggmap(map) + ### Plot 2
  geom_raster(data = dataset, aes(x = Var1, y = Var2, fill = z), alpha = 0.5,  interpolate = TRUE) +
  scale_fill_viridis_c(option = "magma", na.value = "transparent") +
  coord_equal()

The first map looks perfect: neat, clean, lines are defined. First Map

The second one, having added the geom_raster layer, looks (besides wider) slightly blurred, not that crisp. See that the geom_raster line adds a whitish layer ontop of the map (if you look closely it does not even cover it totally). It is absolutely awful and I would like to remove it, or, in other words, I would like it to take a "transparent" colour when the tile produced by geom_raster takes a NA value. Second Map

Any ideas?


Solution

  • If I understood correctly the issue, it seems like the NA in the raster are not completely transparent. See if in scale_fill_viridis_c changing to na.value = NA does what you're looking for.

    library(ggmap)
    library(osmdata)
    box <- c(left = 2.075, bottom = 41.325, right = 2.25, top = 41.47)
    map <- get_stamenmap(bbox = box, maptype = "terrain-lines", zoom = 13)
    
    lon_grid <- seq(2.075, 2.25, length.out = 30)
    lat_grid <- seq(41.325, 41.47, length.out = 30)
    grid <- expand.grid(lon_grid, lat_grid)
    z <- c(rep(NA, 30^2/2), rnorm(30^2/2))
     <- cbind(grid, z)
    
    ggmap(map) + ### Plot 2
      geom_raster(data = dataset, aes(x = Var1, y = Var2, fill = z), alpha = 0.5,  interpolate = TRUE) +
      scale_fill_viridis_c(option = "magma", na.value = NA) +
      coord_equal()
    
    

    This is what it outputs: enter image description here