ggplot2 map where shape and color represent two groups

I have a dataset of work location (i.e., latitude and longitude), if the study system is freshwater, saltwater (i.e., marine), or both, and the home country of the investigator. I am trying to display this information on a world map such that the shape of each point represents the system (i.e., System) and the color fill of each shape represents home country (i.e., Home). I can get this information plotted correctly on the map however, the legend for home country is not correct. How can I make the legend points for home country represent the assigned color?

library("ggplot2")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

rm(list = ls())
dev.off()

# Example dataset
people2 <- as.data.frame(matrix(ncol = 4, nrow = 3))
colnames(people2)[1:4] <- c('Lat','Lon','System','Home')
people2$Lat <- c(46.86197, 45.83034, -27.96287)
people2$Lon <- c(-113.98460, 15.98697, 153.38142)
people2$System <- c('Both', 'Freshwater', 'Marine')
people2$Home <- c('China', 'Croatia', 'Australia')
people2 <- st_as_sf(people2, 
                   coords = c("Lon", "Lat"),
                   crs = 4326)

# Geographic representation of Earth's countries
world <- ne_countries(scale = "medium", returnclass = "sf")

# Geographic representation of Earth's large lakes
lakes <- rnaturalearth::ne_download(scale = 110,
                                    type = 'lakes',
                                    category = 'physical') %>%
  sf::st_as_sf(lakes110, crs = epsg:4326)

# Geographic representation of Earth's oceans
ocean <- st_polygon(list(cbind(c(seq(-180, 179, len = 100), rep(180, 100), 
                                 seq(179, -180, len = 100), rep(-180, 100)),
                               c(rep(-90, 100), seq(-89, 89, len = 100),
                                 rep(90, 100), seq(89, -90, len = 100))))) |>
  st_sfc(crs = "WGS84") |>
  st_as_sf()

# Global basemap in Robinson projection
p0 <- ggplot(data = world) +
  geom_sf(data = ocean, fill = "aliceblue") +
  geom_sf(color = "gray60", fill = "white") +
  geom_sf(data = lakes,
          mapping = aes(geometry = geometry),
          color = "gray60",
          fill = "aliceblue") +
  coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", expand = F) +
  theme_bw() +
  theme(panel.border = element_blank(),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())
p0

# Colors, shapes, and shape sizes for map plotting
cols <- c("Australia" = "#a50026",
           "China" = "#fee090",
           "Croatia" = "#4575b4")
shapes <- c("Freshwater" = 21, "Marine" = 22, "Both" = 24)
shape_size = c("Freshwater" = 2.5, "Marine" = 2.2, "Both" = 2.2)

# Plotting the geographic data
p0 +
  geom_sf(data = people2,
          aes(fill = Home,
              shape = System,
              size = System)) +
  scale_fill_manual(values = cols) + 
  scale_shape_manual(values = shapes) +
  scale_size_manual(values = shape_size) +
  coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", expand = F) +
  guides(fill = guide_legend(ncol = 1, title.position = 'top', title.hjust = 0.5),
         shape = guide_legend(ncol = 1, title.position = 'top', title.hjust = 0.5)) +
  labs(fill = "Home Country") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = 'bottom',
        legend.title = element_text(size = 14, color = 'black'),
        legend.text = element_text(size = 14),
        legend.background = element_blank())

NOTE how the location shape and color are displayed correctly on the map BUT the home country legend shows black points rather than color-filled points. How can I make the home country legend show filled color points?

>Solution :

You can set the shape that will be used for fill aesthetic in the legend to be a coloured one using (eg):

guides(fill = guide_legend(override.aes=list(shape = 21)))

Full code for image:

# Example dataset
people2 <- as.data.frame(matrix(ncol = 4, nrow = 3))
colnames(people2)[1:4] <- c('Lat','Lon','System','Home')
people2$Lat <- c(46.86197, 45.83034, -27.96287)
people2$Lon <- c(-113.98460, 15.98697, 153.38142)
people2$System <- c('Both', 'Freshwater', 'Marine')
people2$Home <- c('China', 'Croatia', 'Australia')
shapes <- c("Freshwater" = 21, "Marine" = 22, "Both" = 24)
ggplot(data=people2 , aes(x=Lat,y=Lon, shape=System, fill=Home )) + 
  geom_point() + 
  scale_shape_manual(values = shapes) + 
  guides(fill = guide_legend(override.aes=list(shape = 21)))