R: how to create a heat map of averaged values from a grid and plot it with ggplot?

Question

I have a data frame (see below) with over 50 000 values, each associated to a position (lat, lon). I would like to calculate the average value for each cell of a 5° latitude x 5° longitude grid in order to create a heat map. The final goal is to plot this grid over a bathymetry map.

I looked at similar questions like this one Average values of a point dataset to a grid dataset. But I couldn't replicate these examples with my own data. Saddly, I am stuck at the first step which is creating the grid.

My data look like this:

library(sp)
library(proj4)

coordinates(data) <- c("lon", "lat")        
proj4string(data) <- CRS("+init=epsg:4326") #defined CRS to WGS 84
df<- data.frame(data)

> head(df)
         lon      lat  value
1 -48.1673562 57.71791  822.9
2 -48.7430053 57.83568 1302.3
3 -48.5662663 57.82087 1508.0
4 -48.3252052 58.29815  224.0
5 -47.1716772 58.42417   38.0
6 -46.4098311 58.67651  431.2
7 -45.8071218 58.70022  365.6
8 -45.5558936 58.46975   50.0

Ideally, I would like to plot the grid on a map from the marmap package using ggplot2 (see below):

library(marmap)
library(ggplot2)

atlantic <- getNOAA.bathy(-80, 40, 0, 90, resolution = 25, keep = TRUE)

atl.df <- fortify(atlantic)

map <- ggplot(atl.df, aes(x=x, y=y)) +
  geom_raster(aes(fill=z), data=atl.df) +
  geom_contour(aes(z=z),
               breaks=0, #contour for continent
               colour="black", size=1) +
  scale_fill_gradientn(values = scales::rescale(c(-5000, 0, 1, 2400)),
                       colors = c("steelblue4", "#C7E0FF", "gray40", "white"))

score 1 · Accepted Answer · answered Sep 10 '20 at 03:20

It sounds like you want to cut your numerical variables (lat & lon) into even intervals and summarise the values within each interval. Does the following work for you?

library(dplyr)

df2 <- df %>%
  mutate(lon.group = cut(lon, breaks = seq(floor(min(df$lon)), ceiling(max(df$lon)), by = 5),
                         labels = seq(floor(min(df$lon)) + 2.5, ceiling(max(df$lon)), by = 5)),
         lat.group = cut(lat, breaks = seq(floor(min(df$lat)), ceiling(max(df$lat)), by = 5),
                         labels = seq(floor(min(df$lat)) + 2.5, ceiling(max(df$lat)), by = 5))) %>%
  group_by(lon.group, lat.group) %>%
  summarise(value = mean(value), .groups = "drop") %>%
  mutate(across(where(is.factor), ~as.numeric(as.character(.x))))

Sample data:

set.seed(444)

n <- 10000
df <- data.frame(lon = runif(n, min = -100, max = -50),
                 lat = runif(n, min = 30, max = 80),
                 value = runif(n, min = 0, max = 1000))

> summary(df)
      lon              lat            value          
 Min.   :-99.99   Min.   :30.00   Min.   :   0.1136  
 1st Qu.:-87.55   1st Qu.:42.45   1st Qu.: 247.2377  
 Median :-75.29   Median :55.11   Median : 501.4165  
 Mean   :-75.12   Mean   :55.01   Mean   : 499.5385  
 3rd Qu.:-62.69   3rd Qu.:67.63   3rd Qu.: 748.8834  
 Max.   :-50.01   Max.   :80.00   Max.   : 999.9600

Comparison of before & after data:

gridExtra::grid.arrange(
  ggplot(df, 
         aes(x = lon, y = lat, colour = value)) + 
    geom_point() + 
    ggtitle("Original points"),
  ggplot(df2, 
         aes(x = lon.group, y = lat.group, fill = value)) + 
    geom_raster() + 
    ggtitle("Summarised grid"),
  nrow = 1
)

Thank you so much @Z.Lin for your answer. This is precisely what I was looking for as it allows to specifiy the size of the grid cells directly. — Don-Jean, Sep 25 '20 at 15:03

score 0 · Answer 2 · answered Sep 10 '20 at 06:38

As (almost!) always, there's a function for that. I believe marmap::griddify() is what you are looking for. The help file states:

Transforms irregularly spaced xyz data into a raster object suitable to create a bathy object with regularly spaced longitudes and latitudes.

Here's a script using your coordinates:

library(marmap)
library(ggplot2)

# Create fake data
set.seed(42)
n <- 10000
data_irregular <- data.frame(lon = runif(n, min = -80, max = 40),
                             lat = runif(n, min = 0, max = 90),
                             value = runif(n, min = 0, max = 1000))

# Fit data into a grid of 30 cells in longitude and 50 cells in latitude
data_grid <- as.bathy(griddify(data_irregular, nlon = 30, nlat = 50))
fortified_grid <- fortify(data_grid)

# Get bathymetric data to plot continent contours
atlantic <- getNOAA.bathy(-80, 40, 0, 90, resolution = 25)
atl_df <- fortify(atlantic)

# Plot with ggplot with gridded data as tiles
map <- ggplot(atl_df, aes(x = x, y = y)) +
  geom_raster(data = fortified_grid, aes(fill = z)) +
  geom_contour(data = atl_df, aes(z = z), 
               breaks = 0, # contour for continent
               colour = "black", size = 1) +
  scale_fill_gradientn(values = scales::rescale(c(-5000, 0, 1, 2400)),
                       colors = c("steelblue4", "#C7E0FF", "gray40", "white")) +
  labs(x = "Longitude", y = "Latitude", fill = "Value")


map +
  theme_bw()

And here is the result:

Thank you @Benoit. Your answers are always very helpful. This works well and the details on how to integrate the results in ggplot2 were very useful. — Don-Jean, Sep 25 '20 at 15:05

R: how to create a heat map of averaged values from a grid and plot it with ggplot?

2 Answers2