Search code examples
rheatmapspatialrgdal

Spatial heatmap with given value for colour

I built upon this thread to plot a spatial heatmap. The difference is that I do not want to compute the density of points because I already have the value for the level of "heat". In detail I want to plot the population density of the canton of Berne (Switzerland) with colour gradients.

The population data is from the Swiss statistical office and counts the inhabitants per hectare (100m x 100m squares), downloadable here (the "STATPOP2020.csv" file). Building on jlhoward answer in thread my code so far is:

library(tidyverse)
library(rgdal)
library(GADMTools)
library(RColorBrewer) 

# conversion from swiss coordinate system to wgs-84
lv03_wgs_lat <- function (y, x){
  y_aux <- (y - 600000)/1000000
  x_aux <- (x - 200000)/1000000
  lat <- {16.9023892 +
      3.238272 * x_aux -
      0.270978 * (y_aux^2) -
      0.002528 * (x_aux^2) -
      0.0447   * (y_aux^2) * x_aux -
      0.0140   * (x_aux^3)}
  lat <- lat * 100/36
  return(lat)  
}

lv03_wgs_lon <- function (y, x){
  y_aux <- (y - 600000)/1000000
  x_aux <- (x - 200000)/1000000
  lon <- {2.6779094 +
      4.728982 * y_aux +
      0.791484 * y_aux * x_aux +
      0.1306   * y_aux * (x_aux^2) -
      0.0436   * (y_aux^3)}
  lon <- lon * 100/36
  return(lon)
}

# read in data
d_pop <- read_csv2("STATPOP2020.csv") %>%
  select(1:6) %>%
  rename(TOT = 6) %>%
  mutate(lon = lv03_wgs_lon(X_KOORD, Y_KOORD),
         lat = lv03_wgs_lat(X_KOORD, Y_KOORD))

# filter swiss data to canton of berne
d_map_ch <- gadm_sf.loadCountries("CHE", level = 1)
d_map_be <- gadm_subset(d_map_ch, level = 1, regions = "Bern", usevar = NULL)
d_map_points <- st_as_sf(d_pop, coords = c("lon", "lat"), crs = 4326)
d_pop_be <- bind_cols(d_pop,
                      as_tibble(t(st_contains(x = d_map_be$sf, y = d_map_points, sparse = FALSE))) %>%
                        rename(in_be = V1)) %>%
  filter(in_be)

# building on previous answer, administrative border is disregarded
ggplot(d_pop_be, aes(x = E_KOORD, y = N_KOORD)) + 
  stat_density2d(aes(fill = TOT), alpha = 0.4, geom = "polygon")+
  scale_fill_gradientn(colours=rev(brewer.pal(7,"Spectral")))+
  xlim(2555000, 2678000) +
  ylim(1130000, 1245000) +
  coord_fixed()

As the heatmap gets plotted, I do not get a legend for the "heat", population density, and I am not able to add colours to the gradients. Does anybody know how to add those?

Solution

  • The problem, as you have already established, is that you want a contour map that represents population density, not the density of measurements, which is what stat_density_2d does. It is possible to create such an object in R, but it is difficult when the measurements are not spaced regularly on a grid (as is the case with this data). It may be best to use geom_point here for that reason:

    ggplot(d_pop_be, aes(x = E_KOORD, y = N_KOORD)) + 
  geom_point(aes(color = log(TOT), alpha = exp(TOT))) +
  scale_colour_gradientn(colours=rev(brewer.pal(7,"Spectral")),
                         breaks = log(c(1, 10, 100, 1000)),
                         labels = c(1, 10, 100, 1000),
                         name = "Population density\n(People per hectare)")+
  xlim(2555000, 2678000) +
  ylim(1130000, 1245000) +
  guides(alpha = guide_none()) +
  coord_fixed()

    enter image description here

    If you want a filled contour you will have to manually create a matrix covering the area of interest, get the mean population in each bin, convert that into a data frame, then use geom_contour_filled:

    z <- tapply(d_pop_be$TOT, list(cut(d_pop_be$E_KOORD, 200), 
                               cut(d_pop_be$N_KOORD, 200)), mean, na.rm = TRUE)

df <- expand.grid(x = seq(min(d_pop_be$E_KOORD), max(d_pop_be$E_KOORD), length = 200),
                  y = seq(min(d_pop_be$N_KOORD), max(d_pop_be$N_KOORD), length = 200))


df$z <- c(tapply(d_pop_be$TOT, list(cut(d_pop_be$E_KOORD, 200), 
                  cut(d_pop_be$N_KOORD, 200)), mean, na.rm = TRUE))

df$z[is.na(df$z)] <- 0

 ggplot(df, aes(x, y)) + 
 geom_contour_filled(aes(z = z), breaks = c(1, 5, 20, 50, 100, 1000)) +
 scale_fill_manual(values = rev(brewer.pal(5, "Spectral")))

    enter image description here