I am trying to overlay lunge_count averages across hexagons on a map; however, I am only able to produce two separate plots (one of the geom_hex and one of the map).

#libraries
library(dplyr) # for tidying
library(lubridate) # for working with datetimes 
library(readr) # for reading in data
library(sf) # for working with spatial data 
library(ggplot2) # for plotting/mapping 
library(ggspatial) # mapping aesthetics 
library(rnaturalearth) # for getting land data
library(rnaturalearthhires) # for getting land data
library(patchwork) #patching maps together
library(leaflet) # Has neat colour palette function
library(plotrix)

Here is the data I am using:

      ptt depth_m duration_s duration_m lunge_count year               id       lon      lat
9946 5883      51        778   12.96667           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9947 5883      53        714   11.90000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9948 5883      49        696   11.60000           1 2023 2023MX-Mno-05883 -107.2284 22.81825
9949 5883      45        756   12.60000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9950 5883      45        780   13.00000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9951 5883      45        752   12.53333           0 2023 2023MX-Mno-05883 -107.2284 22.81825

Here is the code I am using:

#ocean map information

esri_ocean <- paste0('/','Ocean/World_Ocean_Base/MapServer/tile/${z}/${y}/${x}.jpeg')

#create ocean basemap

basemap36 <- ggplot() + annotation_map_tile(type = esri_ocean, zoomin = 1, progress = "none") + coord_sf(crs = 4326,xlim = c(-100, -125),ylim = c(17, 36))

#plotting map and hexagon lunge count averages

basemap36 & ggplot(breeding_redo_5801, aes(lon, lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result

#other attempts

I was able to get the plot to overlay the map with the code below; however, it produced sums instead of averages across the hexagons.

basemap36 & geom_hex(data = breeding_redo_5801, aes(x = lon, y = lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result of this attempt

I am trying to overlay lunge_count averages across hexagons on a map; however, I am only able to produce two separate plots (one of the geom_hex and one of the map).

#libraries
library(dplyr) # for tidying
library(lubridate) # for working with datetimes 
library(readr) # for reading in data
library(sf) # for working with spatial data 
library(ggplot2) # for plotting/mapping 
library(ggspatial) # mapping aesthetics 
library(rnaturalearth) # for getting land data
library(rnaturalearthhires) # for getting land data
library(patchwork) #patching maps together
library(leaflet) # Has neat colour palette function
library(plotrix)

Here is the data I am using:

      ptt depth_m duration_s duration_m lunge_count year               id       lon      lat
9946 5883      51        778   12.96667           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9947 5883      53        714   11.90000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9948 5883      49        696   11.60000           1 2023 2023MX-Mno-05883 -107.2284 22.81825
9949 5883      45        756   12.60000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9950 5883      45        780   13.00000           0 2023 2023MX-Mno-05883 -107.2284 22.81825
9951 5883      45        752   12.53333           0 2023 2023MX-Mno-05883 -107.2284 22.81825

Here is the code I am using:

#ocean map information

esri_ocean <- paste0('https://services.arcgisonline.com/arcgis/rest/services/','Ocean/World_Ocean_Base/MapServer/tile/${z}/${y}/${x}.jpeg')

#create ocean basemap

basemap36 <- ggplot() + annotation_map_tile(type = esri_ocean, zoomin = 1, progress = "none") + coord_sf(crs = 4326,xlim = c(-100, -125),ylim = c(17, 36))

#plotting map and hexagon lunge count averages

basemap36 & ggplot(breeding_redo_5801, aes(lon, lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result

#other attempts

I was able to get the plot to overlay the map with the code below; however, it produced sums instead of averages across the hexagons.

basemap36 & geom_hex(data = breeding_redo_5801, aes(x = lon, y = lat, z = lunge_count)) & stat_summary_hex(fun = mean)

#result of this attempt

• r
• ggplot2
• mapping

2 Answers 2

I believe this will work, but I don't think we have enough example data to try it:

basemap36 & stat_summary_hex(
  data = breeding_redo_5801, 
  mapping = aes(lon, lat, z = lunge_count), 
  fun = mean
)

I think the mistake is here:

basemap36 & ggplot(breeding_redo_5801, aes(lon, lat, z = lunge_count)) & stat_summary_hex(fun = mean)

You should use the operator + to add layers to a ggplot object. See a minimal reprex where I simulate breeding_redo_5801:

# Minimal libraries to use
library(ggplot2)
library(ggspatial)
library(dplyr)

esri_ocean <- paste0('https://services.arcgisonline.com/arcgis/rest/services/','Ocean/World_Ocean_Base/MapServer/tile/${z}/${y}/${x}.jpeg')


basemap36 <- ggplot() + annotation_map_tile(type = esri_ocean, zoomin = 1, progress = "none") + coord_sf(crs = 4326,xlim = c(-100, -125),ylim = c(17, 36))

Now simulate something like (minimal) breeding_redo_5801 to get more points:

breeding_redo_5801 <- tibble::tibble(
  lon = runif(10000, min = -110, max = -105),
  lat = runif(10000, min = 20, max = 25),
  lunge_count = runif(10000, min = 1, max = 5)
) %>%
  arrange(lon, lat)

breeding_redo_5801
#> # A tibble: 10,000 × 3
#>      lon   lat lunge_count
#>    <dbl> <dbl>       <dbl>
#>  1 -110.  20.1        2.40
#>  2 -110.  21.5        2.60
#>  3 -110.  22.6        2.88
#>  4 -110.  24.5        1.71
#>  5 -110.  24.3        4.32
#>  6 -110.  21.7        4.44
#>  7 -110.  22.6        2.16
#>  8 -110.  20.6        4.67
#>  9 -110.  24.9        3.82
#> 10 -110.  22.7        3.91
#> # ℹ 9,990 more rows

And by last add your layer:

basemap36 +
  # Add hex layer
  stat_summary_hex(data = breeding_redo_5801, aes(lon, lat, z = lunge_count), 
                   fun = mean)