How can I get multiple icon types on ggplot?

Question

Hopefully this is a very simple solution. I'm completely new to R Studio and don't even completely understand what code I DO have so far. So please let me know if I need to give more info.

I'm mapping locational data for different plant species, coloring them by decade. I also need to put triangles for plants in an herbarium and circles for the ones observed in the wild. This is what I have so far, along with the image it produces.

In my data, I have a column titled "basisOfRecord" which has the two categories of either human observation or preserved specimen. How would I change my code to assign the map icons? I kind of had this project thrown at me and have been messily googling without a foundation in R or even coding in general, so I'm sorry if my script is really messy.

Here is a link to the map data - https://www.dropbox.com/scl/fo/v2v8x3zy4uhlw3mtan28l/h?rlkey=60g8s28z1mix9co2podjkb16p&dl=0

And here is a link to the CSV - https://www.dropbox.com/scl/fi/9t7qpxo1113uban1qpsib/americanus.csv?rlkey=nguai55lmlsc23i48h1j2zeee&dl=0

library(tidyverse)
library(lubridate)
library(viridis)
library(sf) 


americanus <- read.csv("C:/Users/MayeJ/Documents/Map Data/GBIF/Americanus/americanus.csv")  

month <- 1  
day <- 1
americanus$TrueDate <- as.Date(paste(americanus$TrueDate, month, day, sep = "-"))
print(americanus$TrueDate)

americanus <- americanus[complete.cases(americanus$decimalLatitude, americanus$decimalLongitude), ]

obs.year = lubridate::year(americanus$TrueDate) 
decade.cha = as.character(floor(obs.year / 10) * 10) 
decade = (floor(obs.year / 10) * 10)
americanus <- st_as_sf(americanus, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

crs = "WGS84"
agr = "constant"
remove = F



class(americanus)


north.america.bound <- st_read("/Users/MayeJ/Documents/Map Data/NA_PoliticalDivisions/NA_PoliticalDivisions/data/bound_l/boundary_l_v2.shp")



americanus.map <- ggplot(data = north.america.bound) + 

  theme_bw() + 
  geom_sf(fill = "gray90", color = "black") + 
  geom_sf(data = americanus, 
          size = 3, 
          shape = 21, 
          alpha = 0.7, 
          aes(fill = decade.cha)) + 
  scale_fill_viridis(option = "B", discrete = T, direction = 1) +

  coord_sf(xlim = c(-3235000, 3062000),
           ylim = c(-900000, 4007000)) 


americanus.map



ggsave(plot = americanus.map,
       "C:/Users/MayeJ/Documents/Map Data/GBIF/Americanus/americanusmap.jpg",
       device = "jpg", dpi = 300, units = "in",
       width = 10, height = 7.5)

Here is the map

I've tried looking for examples of maps that use ggplots and sf, to see if they use more than one symbol at a time, but haven't found anything yet. I just don't have any kind of formatting intuition with coding, so I wouldn't even know where to start in experimenting.

Answer 1

The issue is that your basisOfRecord column has six categories. Four are simply different versions of PreservedSpecimen (e.g. there is also "preservedspecimen"), then there is a "HumanObservation" and the last one is called "voucher". Hence, when providing only two values for shape (as I suggested in the comment) you will get an error while when not including scale_shape_manual you get the first six default shapes which support only the color aes but not fill.

In the code below I recoded the four variants of "PreservedSpecimen" as one category "PreservedSpecimen" and added another shape value for "voucher". However, there was still an issue with fill colors not showing up in the legend which I fixed using guides(fill = guide_legend(override.aes = list(shape = 21))). Additionally I simplified your code a bit, especially the computation of the decades and finally I put the legend at the bottom.

library(tidyverse)
library(sf)

americanus <- read.csv("americanus.csv")
americanus <- americanus[c("TrueDate","decimalLatitude", "decimalLongitude", "basisOfRecord")]

americanus$decade.cha <- factor(10 * americanus$TrueDate %/% 10)

americanus <- americanus[complete.cases(americanus$decimalLatitude, americanus$decimalLongitude), ]

americanus <- st_as_sf(americanus, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

north.america.bound <- st_read("PoliticalBoundaries_Shapefiles/NA_PoliticalDivisions/data/bound_l/boundary_l_v2.shp")

americanus$basisOfRecord[grepl("^(P|p)reserve", americanus$basisOfRecord)] <- "PreservedSpecimen"
ggplot(north.america.bound) +
  geom_sf(fill = "gray90", color = "black") +
  geom_sf(
    data = americanus,
    aes(fill = decade.cha, shape = basisOfRecord),
    size = 3,
    alpha = 0.7
  ) +
  scale_shape_manual(values = c(21, 24, 22)) +
  scale_fill_viridis_d(option = "B", direction = 1) +
  coord_sf(
    xlim = c(-3235000, 3062000),
    ylim = c(-900000, 4007000)
  ) +
  theme_bw() +
  theme(legend.position = "bottom", legend.box = "vertical") +
  guides(fill = guide_legend(override.aes = list(shape = 21)))