In the last section we discussed how to obtain and then process data from the Global Biodiversity Information Facility using the rgbif package from rOpenSci in R with RStudio.

In this section we explore one of the real strengths of GBIF data by visualizing the data for Kenya on a map using the popular leaflet library. We will then create an online interactive version of the map as a Shiny app.

The image below from the GBIF website displays a global map of occurrence records.

It is important to bear in mind that we can only visualize those records for a particular country or region where georeferenced coordinates (latitude and longitude) are provided with GBIF occurrence records. As such records with coordinates will be a subset of the wider occurrence records set. In connection with access and benefit-sharing under the Nagoya Protocol this is one, among other, reasons why researchers should be encouraged to record georeferenced data when collecting samples.

Most visual representations are static, with a growing trend towards maps that we can zoom into (as in the case of GBIF maps). However, this raises the question of whether it is possible to create interactive maps that allow the reader to explore the data in more detail.

A growing number of software tools are now available for the creation of interactive maps. In this section we will use the popular open source leaflet JavaScript library to create interactive maps for species data that are interactive and mobile friendly.

The creation of maps with unfamiliar software can be somewhat intimidating. We will use the leaflet package developed by RStudio and use the RStudio tutorial as a guide to using leaflet and adapt as we go along.

Happily, it is quite straightforward to learn to use leaflet and the data preparation steps in this walk through could quite easily be done in either another programming language or with Open Office Calc or Excel. Creating an online app with Shiny is more challenging but is easier when the basic principles behind Shiny are understood.

Getting Started

First, make sure that you have R and RStudio installed for your system by following the installation instructions.

Next, install the leaflet package and if you have not done so already install the tidyverse packages. tidyverse packages are a set of packages that focus on data wrangling and are extremely useful for data preparation.


Now load the libraries we will need. Don’t worry if you see a warning message.


Getting the GBIF Data

Let’s try to visualize some GBIF data occurrence data. For this test we will use Kenya occurrence data from January 2017 with the following doi for citation: doi:10.15468/dl.b04fyt.

An earlier article walked through the process of downloading an occurrence dataset from GBIF and then cleaning up the data to remove those records that do not contain coordinate data. You can find that article here.

The file containing those results can be found on Github here.


Preparing the dataset for mapping

There are a number of steps that we need to take to prepare the dataset for mapping. These involve:

  1. Renaming columns
  2. Filtering the data for species records (see the previous gbif article)
  3. Ensuring that the georeferenced coordinate columns are complete (see the previous gbif article)

Renaming columns

The leaflet package will attempt to guess the name of the latitude and longitude columns in a dataset. Because GBIF calls them decimallatitude and decimallongitude, which is difficult to guess, it makes sense to rename them using the nice and easy dplyr rename() function.

kenya_occurrence <- dplyr::rename(kenya_occurrence, latitude = decimallatitude, 
    longitude = decimallongitude)

Filtering the data

In the previous section we filtered the data on the taxonrank column for species. In exploring the data we discovered that the occurrence records are heavily dominated by observations for birds as we can see below.

kenya_gbif %>% drop_na(class) %>% count(class, sort = TRUE)
## # A tibble: 122 × 2
##                 class      n
##                 <chr>  <int>
## 1                Aves 416835
## 2             Insecta  68932
## 3            Mammalia  64124
## 4       Magnoliopsida  48989
## 5          Liliopsida  23984
## 6            Amphibia  14510
## 7            Reptilia  13383
## 8      Actinopterygii   8964
## 9          Gastropoda   4812
## 10 Protosteliomycetes   3939
## # ... with 112 more rows

We now want to filter out these records because they will flood the map. Here we will filter the kenya_occurrence dataset. We will save this dataset to a new set called kenya_occ

kenya_occ <- filter(kenya_occurrence, class != "Aves")

This reduces our occurrence data to 105,853 occurrence records from 444,228.

Creating a map with leaflet

Next, let’s simply adapt the code from the RStudio tutorial to see what we learn and adapt as we go along. This example uses pipes ( %>% for this then that) so that we only have to name the dataset in the first line. If you are new to pipes then try Sean Anderson’s tutorial or the excellent new R For Data Science by Garrett Grolemund and Hadley Wickham.

We will start with a smaller dataset with 20 rows to learn with and then work our way up.

kenya_20 <- kenya_occ[1:20, ]

To get started we use the leaflet() function and run a test with the kenya_20 dataset. Then we will add the default OpenStreetMap tiles used in leaflet. Next we will use addMarkers() and specify longitude and latitude with a pop up marker by selecting the kenya_occurrence species column using the $ dollar sign.

This provides us with some important clues on what will happen when we adapt the RStudio code.

k <- leaflet::leaflet(kenya_20) %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(~longitude, ~latitude, popup = kenya_20$species)

Hmmm, we now have a map with a few points but not the 20 we were expecting. If we click on a blue marker the species name will pop up. So, the code is picking up some of the names but not others. This suggests we need to look at the data in more detail.

If you are working with your own data and have not previously prepared the data by filtering on taxonrank species and to filter NAs, part of the explanation will lie with the need for data preparation. However, it also becomes clear that for some coordinates we have exactly the same georeferences as we see below for Parmotrema andinum.

kenya_20 %>% select(species, latitude, longitude) %>% print()
## # A tibble: 20 × 3
##                      species latitude longitude
##                        <chr>    <dbl>     <dbl>
## 1         Bulbothrix isidiza   -1.274        37
## 2    Caloplaca ochraceofulva   -0.356        36
## 3   Caloplaca flavorubescens   -3.367        39
## 4  Canoparmelia nairobiensis   -1.274        37
## 5        Canoparmelia texana   -0.276        36
## 6     Cerastium afromontanum   -0.164        37
## 7     Cerastium afromontanum   -0.628        37
## 8         Cladonia hedbergii   -0.167        37
## 9         Crassula granvikii   -0.182        37
## 10        Lobelia aberdarica   -0.515        35
## 11        Lobelia gregoriana   -0.628        37
## 12        Lobelia gregoriana   -0.350        37
## 13  Parmelinopsis horrescens   -0.366        35
## 14        Parmotrema andinum   -1.274        37
## 15        Parmotrema andinum   -1.274        37
## 16        Parmotrema cooperi   -0.163        36
## 17        Parmotrema eunetum    1.100        37
## 18      Parmotrema natalense   -0.083        35
## 19      Parmotrema natalense   -0.083        35
## 20      Rinodina groenhartii   -0.417        36

We do not need to take any action in these cases and leaflet is behaving as expected. Let’s sample 500 records across the table and see how it goes.

kenya_500 <- kenya_occ %>% sample_n(500)
## # A tibble: 500 × 44
##        gbifid                           datasetkey
##         <int>                                <chr>
## 1   183706538 09baff40-d4a3-11dc-8719-b8a03c50a862
## 2    47414991 85818aea-f762-11e1-a439-00145eb45e9a
## 3  1270009958 7e380070-f762-11e1-a439-00145eb45e9a
## 4   177626402 68efc55e-f762-11e1-a439-00145eb45e9a
## 5  1056907247 7e380070-f762-11e1-a439-00145eb45e9a
## 6  1052778786 e6c97f6e-e952-11e2-961f-00145eb45e9a
## 7   142746983 b276cb50-d3ea-11dc-ab69-b8a03c50a862
## 8  1230338192 1aaec653-c71c-4695-9b6e-0e26214dd817
## 9   476568766 4bfac3ea-8763-4f4b-a71a-76a6f5f243d3
## 10  912192321 cd6e21c8-9e8a-493a-8a76-fbf7862069e5
## # ... with 490 more rows, and 42 more variables: occurrenceid <chr>,
## #   kingdom <chr>, phylum <chr>, class <chr>, order <chr>, family <chr>,
## #   genus <chr>, species <chr>, infraspecificepithet <chr>,
## #   taxonrank <chr>, scientificname <chr>, countrycode <chr>,
## #   locality <chr>, publishingorgkey <chr>, latitude <dbl>,
## #   longitude <dbl>, coordinateuncertaintyinmeters <dbl>,
## #   coordinateprecision <dbl>, elevation <dbl>, elevationaccuracy <dbl>,
## #   depth <dbl>, depthaccuracy <dbl>, eventdate <chr>, day <int>,
## #   month <int>, year <int>, taxonkey <int>, specieskey <int>,
## #   basisofrecord <chr>, institutioncode <chr>, collectioncode <chr>,
## #   catalognumber <chr>, recordnumber <chr>, identifiedby <chr>,
## #   license <chr>, rightsholder <chr>, recordedby <chr>, typestatus <chr>,
## #   establishmentmeans <chr>, lastinterpreted <chr>, mediatype <chr>,
## #   issue <chr>

We now input kenya_500 into the code.

species <- leaflet(kenya_500) %>% addTiles() %>% addMarkers(~longitude, ~latitude, 
    popup = kenya_500$species)

Our kenya_500 dataset reveals that some of the coordinates are outside the boundaries of Kenya. Let’s take a look at a couple of those records that illustrate the nature of the problem.

kenya_500 %>% filter(species == "Medicago laciniata" | species == "Chaetodon guttatissimus") %>% 
    select(species, latitude, longitude, countrycode, locality, issue) %>% print()
## # A tibble: 2 × 6
##                   species latitude longitude countrycode
##                     <chr>    <dbl>     <dbl>       <chr>
## 1 Chaetodon guttatissimus      -40        40          KE
## 2      Medicago laciniata        0         0          KE
## # ... with 2 more variables: locality <chr>, issue <chr>

We can see that in practice one record record actually lacks coordinates (recorded as 0) while in the other case, the latitude appears to be incorrect. GBIF notes these potential problems in the issue column for country coordinate mismatch and zero coordinate. As such we have a data error.

kenya_500 %>% filter(species == "Medicago laciniata" | species == "Chaetodon guttatissimus") %>% 
    select(issue) %>% print()
## # A tibble: 2 × 1
##                                                                         issue
##                                                                         <chr>

There are a number of ways we could deal with this. We could filter the results with no actual coordinates out and then seek to identify relevant Geodetic data and coordinate issues in the issue column. This would involve unpacking that column (using dplyr::separate_rows() with ; as the separator), then carefully assessing the entries and filtering the results accordingly (e.g. to take out zero_coordinate and country_coordinate_mismatch).

An alternative approach is to adjust the boundaries of the map to exclude records outside Kenya. We will demonstrate this approach here because it is useful for narrowing down maps for other purposes.

Narrowing the Map using boundaries

To narrow the occurrence records to those within the boundaries of Kenya we first need to obtain georeference information on the boundaries. To do that we can use the geonames web service through the geonames R package from Barry Rowlingson.


and then load the library.

## No geonamesUsername set. See and set one with options(geonamesUsername="foo") for some services to work

To use the API we need to sign up for a free account with geonames from here. As we can see above, when we load the geonames package we will receive a warning that no geonamesUsername set. To do that we run this code with your own username.

options(geonamesUsername = "yourusername")

Next, and optionally, run the test files below (generating a lot of returns and not shown) to check that you can access the service.

source(system.file("tests", "testing.R", package = "geonames"), echo = TRUE)

To get the boundaries for a country we need to access the country info. The data that we are after is on this page.

To access that information with the geonames package we use the two letter country code from here or, if you are planning to do a lot of work with country codes, see the useful countrycode package:

kenya_geonames_info <- GNcountryInfo(country = "KE")

The information that we want is the following:

bounds <- kenya_geonames_info %>% select(south, north, east, west)
##       south    north      east      west
## 1 -4.678047 5.019938 41.899078 33.908859

The raw xml from the country table can be seen below.

Next we want to set the boundaries. RStudio helpfully points out how to do this with the following code.

First we can use setView to centre the View using Nairobi as the example.

nairobi <- leaflet() %>% addTiles() %>% setView(lng = 36.821946, lat = -1.292066, 
    zoom = 12)

Note above that longitude (lng) is followed by latitude (lat). If you find yourself unexpectedly in the middle of an ocean then you need to reverse the coordinates.

A useful quick tool for finding the latitude and longitude of a specific place is Otherwise geonames is probably best.

Next, we want to fit the bounds of the map by adding additional information. It is important to note the order in which coordinates are set here (see setView in the leaflet documentation).

The Rstudio examples are:

fitBounds(map, lng1, lat1, lng2, lat2)
# or
setMaxBounds(map, lng1, lat1, lng2, lat2)

Let’s call up the bounds table to see the coordinates.

##       south    north      east      west
## 1 -4.678047 5.019938 41.899078 33.908859

West and East are the longitude coordinates and North and South are the latitudes. Let’s try that with our Nairobi map above.

fb <- fitBounds(nairobi, lng1 = bounds$west, lat1 = bounds$north, lng2 = bounds$east, 
    lat2 = bounds$south)

That works just fine. Now let’s try the setMaxBounds controls.


mb <- setMaxBounds(nairobi, lng1 = bounds$west, lat1 = bounds$north, lng2 = bounds$east, 
    lat2 = bounds$south)

That basically brings us back to the Nairobi view.

If we wanted to clear the boundaries from a particular view we would use clearBounds() to take us back to the world view.

nairobi %>% clearBounds()

Next we would like to see our data inside the fit bounds. First we need our map object to fit the bounds too with the markers and then we need to fit the bounds using fitBounds(). We will use the kenya_500 set.

species <- leaflet(kenya_500) %>% addTiles() %>% addMarkers(~longitude, ~latitude, 
    popup = kenya_500$species) %>% fitBounds(lng1 = bounds$west, lat1 = bounds$north, 
    lng2 = bounds$east, lat2 = bounds$south)

In practice, if we were to zoom out of this map we would still see the other data points outside the boundaries of Kenya. That is fine if we are happy with that.

What we need to do if we want to only see records that fall inside Kenya is to deliberately exclude records that fall outside the boundaries of Kenya.

Excluding coordinates

To exclude coordinates we need to generate information on the max and min latitude and longitude in the kenya_occ_ref dataset. We need four max and min values for that.

max_lat <- max(kenya_occ$latitude)
max_lon <- max(kenya_occ$longitude)
min_lat <- min(kenya_occ$latitude)
min_lon <- min(kenya_occ$longitude)
max_min_lat_lon <- tibble::tibble(max_lat, max_lon, min_lat, min_lon)
## # A tibble: 1 × 4
##   max_lat max_lon min_lat min_lon
##     <dbl>   <dbl>   <dbl>   <dbl>
## 1      62     143     -41    -153

We now need to identify the data points that occur between these coordinates. We will use the dplyr::between() function for that.

lat_between <- kenya_occ$latitude[between(kenya_occ$latitude, -4.678047, 5.019938)]
lon_between <- kenya_occ$longitude[between(kenya_occ$longitude, 33.908859, 41.899078)]

If we print the length of lat_between it is 105552. We can quickly calculate the latitude markers outside the latitude as follows.

length(kenya_occ$latitude) - length(lat_between)
## [1] 301

and the same for longitude

length(kenya_occ$longitude) - length(lon_between)
## [1] 2374

We now want to exclude the lat and long values outside the boundaries. The question then is how to approach this? It makes sense to start by excluding those values that are outside the longitude for Kenya (east, west). Let’s try that using %in% to add a logical TRUE/FALSE where TRUE will be records falling between the longitude values and false will be outside. Then the same for latitude.

kenya_occ$long_out <- kenya_occ$longitude %in% lon_between
kenya_occ$lat_out <- kenya_occ$latitude %in% lat_between

Next we filter to a new data frame on the logical to include first only those where the longitude is within the between values and second to those where the lat is inside the between values. We will overwrite our existing kenya_occ dataset.

kenya_occ <- filter(kenya_occ, long_out == "TRUE") %>% filter(lat_out == "TRUE")

The result is a new data frame with 103295 rows.

Now that we have cleaned up the data let’s create a decent sized sample and then draw a map.

kenya_3000 <- sample_n(kenya_occ, 3000)
map <- leaflet(kenya_3000) %>% addTiles() %>% addMarkers(~longitude, ~latitude, 
    popup = kenya_3000$species) %>% fitBounds(lng1 = bounds$west, lat1 = bounds$north, 
    lng2 = bounds$east, lat2 = bounds$south)

If we zoom out of this map we will see that the vast majority of records falling outside the Kenya boundary box have been excluded. So we now have a map focused on Kenya that excludes extraneous records.

Because we are using a bounding box, and country boundaries are not boxes, we will still encounter extraneous records in nearby countries. To be more precise we would need to find a suitable shapefile for kenya. Or, we could use a geoJSON file such as this one or this one from OpenAfrica (not tested).

For the moment we will test how our new kenya_occ file has worked out and start by reducing the table to relevant columns in a new object we will call kenya_slim.

kenya_slim <- select(kenya_occ, gbifid, kingdom, species, latitude, longitude, 
## # A tibble: 6 × 6
##   gbifid kingdom                   species latitude longitude
##    <int>   <chr>                     <chr>    <dbl>     <dbl>
## 1 653826   Fungi        Bulbothrix isidiza    -1.27        37
## 2 664444   Fungi   Caloplaca ochraceofulva    -0.36        36
## 3 664565   Fungi  Caloplaca flavorubescens    -3.37        39
## 4 673178   Fungi Canoparmelia nairobiensis    -1.27        37
## 5 673247   Fungi       Canoparmelia texana    -0.28        36
## 6 693924 Plantae    Cerastium afromontanum    -0.16        37
## # ... with 1 more variables: locality <chr>

Visualising many thousands of data points is likely to result in memory problems or sluggish performance. We can work our way around this by clustering the markers.


Clustering involves the addition of clusterOptions to the code as below. We will try this with the full kenya_slim dataset.

map <- leaflet(kenya_slim) %>% addTiles() %>% fitBounds(lng1 = bounds$west, 
    lat1 = bounds$north, lng2 = bounds$east, lat2 = bounds$south) %>% addCircleMarkers(~longitude, 
    ~latitude, popup = kenya_slim$species, radius = 1, fillOpacity = 0.5, clusterOptions = markerClusterOptions())

This works fine, and does not create memory problems. If we click on one of the clusters other clusters will open up and we can navigate around the country cluster by cluster. As we drill down we will able to click on the individual markers to see the name of the species.

Changing the Markers

To change the markers from the default we use the addCircleMarkes control along with the specification of the circle radius, weight, opacity, fill and fillOpacity. We will use the kenya_3000 set to illustrate this but you are likely to want to experiment with these settings.

map <- leaflet(kenya_3000) %>% addTiles() %>% addCircleMarkers(~longitude, ~latitude, 
    popup = kenya_3000$species, radius = 1, weight = 2, opacity = 0.5, fill = TRUE, 
    fillOpacity = 0.2)

Using Colours for Categories

In the data above we have changed the data from the large markers to circle markers. Let’s look at this in more detail and use categories to colour the data.

There are a number of controls that can be used to set the colours within leaflet maps including colorBin, colorFactor, colorNumeric and colorQuantile and it is well worth studying the package documentation and the examples in the Leaflet tutorial.

To colour the markers we need to use a colour palette. The ColorBrewer palettes are well known and can be viewed and tested out on the ColorBrewer website. Inside R we can use the RColorBrewer package. First we need to install the package.


Load the library


ColorBrewer contains three types of palette:

  1. Sequential
  2. Diverging
  3. Qualitative

It is the qualitative (categorical) palettes that we want, to give distinct colors to our kingdoms. The qualitative palettes are:

  • Accent 8
  • Dark2 8
  • Paired 12
  • Pastel1 9
  • Pastel2 8
  • Set1 9
  • Set2 8
  • Set3 12

If we visit the color brewer website we can test out the colors across the three types. The name of the palette appears in the URL based on your choices, for example:

We need to define the categories that will be used for the colours. In this case we will want to colour by kingdom. But bear in mind that an individual dataset may not include all of the kingdoms used in the GBIF taxomony. The Kingdoms used in GBIF are:

  • Animalia
  • Archaea
  • Bacteria
  • Chromista
  • Fungi
  • incertae sedis
  • Plantae
  • Protozoa
  • Viruses

So, for use across multiple datasets we will want 9 colour categories. A colour will be allocated to NA (Not Available) by default and so we do not need to worry about that.

With leaflet we can preview the colours we want to apply to our categorical data. This works better run in the console than in an Rmarkdown document.

previewColors(colorFactor("Paired", domain = NULL), LETTERS[1:9])

We now need a vector of names to associate with the colours in the palette. Note that I have cheated by organising the kingdom names to match to the colours that I want. For example, I want Animalia to be red and Plantae to be green (and I am not so bothered about the others).

kingdom <- c("incertae sedis", "Archaea", "Fungi", "Plantae", "Chromista", "Animalia", 
    "Bacteria", "Protozoa", "Viruses")

Then we create a palette from brewer.pal.

my_palette <- brewer.pal(9, "Paired")

Let’s test it.

previewColors(colorFactor(my_palette, levels = kingdom), kingdom)

Now lets create that as a colorFactor.

factpal <- colorFactor(my_palette, levels = kingdom)

Lets add the colours to the maps for our sample data.

species_circle <- leaflet(kenya_3000) %>% addTiles() %>% addCircleMarkers(~longitude, 
    ~latitude, popup = kenya_3000$species, radius = 1, weight = 2, opacity = 0.5, 
    fill = TRUE, fillOpacity = 0.2, color = ~factpal(kingdom))

It will be helpful to add some controls to the map. We can do this by adding layers.

Add Layers

We can readily add layer controls to the map, either to add layers (for example for national parks) or as in this case to allow users to select data by kingdom. The code below is adapted from a Stack Overflow answer from marbel. This approach uses a forloop to loop over the contents of the dataset and assign the data to groups.


groups = unique(kenya_3000$kingdom)

map = leaflet(kenya_3000) %>% addTiles(group = "OpenStreetMap")
for (i in groups) {
    data = kenya_3000[kenya_3000$kingdom == i, ]
    map = map %>% addCircleMarkers(data = data, ~longitude, ~latitude, radius = 1, 
        weight = 2, opacity = 0.5, fill = TRUE, fillOpacity = 0.2, color = ~factpal(kingdom), 
        group = i)
map %>% addLayersControl(overlayGroups = groups, options = layersControlOptions(collapsed = FALSE))

It may be that there is a more efficient (or rather legible) way of going about this but we will stick with this for the time being.

We have now created a control layer. There a quite a few other options we might pursue and we will look at those below. For the moment we will simply look at options to export our map.

Adding Hyperlinked Popups to the Map Points

As we saw earlier, we can simply add popup labels to the data using popup = data$species. However, we can go quite a lot further by creating hyperlinks to make our map interactive. For that we will need a function that will generate the hyperlinks and labels that we might want. We will then need to combine the hyperlinks into a single field for use as the labels.

The function allows for hyperlinks to be constructed for Google, Crossref, GBIF and the Lens (a patent database). The function takes four arguments. The vector of terms to be used in the search, the label for our marker (e.g. Lookup on Google) and the type of link (e.g. Google).

# function to create hyperlinks
map_url <- function(query, label = "NULL", type = "NULL") {
    href <- "<a href="
    close_href <- ">"  #included for flexibility in labelling
    close_a <- "</a>"
    if (type == "google") {
        query <- stringr::str_replace_all(query, " ", "+")
        google_base <- ""
        url <- paste0(google_base, query)
        out <- paste0(href, shQuote(url), close_href, label, close_a)
    if (type == "crossref") {
        query <- stringr::str_replace_all(query, " ", "+%2B")
        crossref_base <- ""
        url <- paste0(crossref_base, query)
        out <- paste0(href, shQuote(url), close_href, label, close_a)
    if (type == "gbif") {
        query <- stringr::str_replace_all(query, " ", "+")
        gbif_base <- ""
        url <- paste0(gbif_base, query)
        out <- paste0(href, shQuote(url), close_href, label, close_a)
    if (type == "lens") {
        # note restriction to main jurisdictions and no stemming to reduce
        # duplication and false positives
        query <- stringr::str_replace_all(query, " ", "+")
        lens_base <- ""
        url <- paste0(lens_base, "%22", query, "%22", "&jo=true&j=EP&j=JP&j=US&j=WO&st=false&n=50")
        out <- paste0(href, shQuote(url), close_href, label, close_a)

Next we need to use the function to add columns to our dataset. Note that here we are using the species name as the query.

# create columns with formatted links
kenya_3000$google <- map_url(kenya_3000$species, label = "Lookup Google", type = "google")
kenya_3000$crossref <- map_url(kenya_3000$species, label = "Lookup Crossref", 
    type = "crossref")
kenya_3000$lens <- map_url(kenya_3000$species, label = "Lookup Patents", type = "lens")
kenya_3000$gbif <- map_url(kenya_3000$species, label = "Lookup GBIF", type = "gbif")

Next we need to combine the columns into one label with the HTML to make them look reasonably neat.

# combine links for use as a popup
kenya_3000$combined_label <- paste0("<br>", "<strong>", kenya_3000$species, 
    "</strong>", "</br>", "<br>", kenya_3000$google, "</br>", "<br>", kenya_3000$gbif, 
    "</br>", "<br>", kenya_3000$crossref, "</br>", "<br>", kenya_3000$lens, 

We now use the combined_label column as the value for the pop up field in our map.


groups = unique(kenya_3000$kingdom)

map = leaflet(kenya_3000) %>% addTiles(group = "OpenStreetMap")
for (i in groups) {
    data = kenya_3000[kenya_3000$kingdom == i, ]
    map = map %>% addCircleMarkers(data = data, ~longitude, ~latitude, popup = data$combined_label, 
        radius = 1, weight = 2, opacity = 0.5, fill = TRUE, fillOpacity = 0.2, 
        color = ~factpal(kingdom), group = i)
map %>% addLayersControl(overlayGroups = groups, options = layersControlOptions(collapsed = FALSE))