SLF occurrence records: GBIF, lydemapR, literature

Overview

In the last vignette, we downloaded bioclimatic rasters from CHELSA and harmonized them to a standard format, including a the resolution, origin and extent.

This vignette will provide instructions to retrieve Lycorma delicatula (SLF) presence records for the future creation of MaxEnt models. MaxEnt is a presence-only modeling software, and so it does not require recorded absence data to predict the suitable area for SLF. Four categories of data sources will be used in this analysis: GBIF (Global Biodiversity Information Facility), lydemapR, various pieces of peer-reviewed literature, and natural history notes. These data will then be tidied, spatially thinned and compiled into a single .rds file of SLF presence records that can be loaded into MaxEnt.

The first step will be to retrieve data from GBIF and lydemapR. GBIF is an open-access platform for biodiversity data that gathers from various databases and citizen science platforms. Data from this source represents globally distributed presences of SLF. The lyde() function within the Lydemapr package gives access to nearly one million SLF presence records within the United States, largely obtained from biological field surveys by various state and federal departments of agriculture.

These data will need to be cleaned and tidied during this step. To get to this tidy dataset, the data will be checked for inconsistent and false records. The data will also be spatially thinned so that points are no less than 10km (the resolution of the climate data, 5 arcminutes) to eliminate the effects of sampling bias.

The second step will be to combine these records with data gathered from peer-reviewed literature and from natural history notes. A majority of these records are from established populations of SLF within China, South Korea and southeast Asia. These records are especially important because China and southeast Asia represent the native range for SLF, and there is very little data on the extent of its native range. These records are also important because MaxEnt correlates presence records with the local climate; thus, if there is very little characterization of the native range for this species, it is unlikely that models for SLF will make realistic predictions for its potential niche elsewhere.

The final step of this analysis will be to organize the different datasets into a single .rds file that can be loaded into MaxEnt. MaxEnt requires that input datasets contain only 3 columns, in this order: species (the scientific name of the species), x (longitude), and y (latitude). Lastly, the data will go through a second round of spatial thinning now that the datasets have been joined.

Setup

First, I will load in the necessary packages.

# general tools
library(tidyverse)  #data manipulation
library(here) #making directory pathways easier on different instances
# here() is set at the root folder of this package
library(devtools) # installing packages not from CRAN

# species occurrence data
library(lydemapr) # field survey data for SLF
library(rgbif) #query gbif and format as a dataframe
library(GeoThinneR) # spatial thinning of species occurrences
library(CoordinateCleaner) # tidying of points

# shapefile data
library(rnaturalearth)
library(rnaturalearthdata)
library(rnaturalearthhires)

# spatial
library(terra)

# aesthetics
library(patchwork) #nice plots
library(knitr) # nice rmd tables

Note: I will be setting the global options of this document so that only certain code chunks are rendered in the final .html file. I will set the eval = FALSE so that none of the code is re-run (preventing files from being overwritten during knitting) and will simply overwrite this in chunks with plots.

I will also download a world map for plotting the data I download from GBIF and lydemapR. I will use the rnaturalearth package to download a world map at a scale of 10m.

# check which types of data are available
# these are in the rnaturalearth package
data(df_layers_cultural) 
# I will use states_provinces

# if the file isnt already downloaded, download it
if (!file.exists(file.path(here::here(), "data-raw", "ne_countries", "ne_10m_admin_0_countries.gpkg"))) {
  
  # get metadata
  ne_metadata <- rnaturalearth::ne_find_vector_data(
    scale = 10,
    category = "cultural",
    getmeta = TRUE
  ) %>%
    dplyr::filter(layer == "admin_0_countries")
  # URL to open metadata
  utils::browseURL(ne_metadata[, 3])
  
  # download
  countries_sf <- rnaturalearth::ne_download(
    scale = 10, # highest resolution
    type = "admin_0_countries", # states and provinces
    category = "cultural",
    destdir = file.path(here::here(), "data-raw", "ne_countries"),
    load = TRUE, # load into environment
    returnclass = "sf" # shapefile
  )
  
  # else, just import
  } else {
    
   countries_sf <- rnaturalearth::ne_load(
      scale = 10,
      type = "admin_0_countries",
      category = "cultural",
      destdir = file.path(here::here(), "data-raw", "ne_countries"),
      returnclass = "sf"
    )
  
  }

Most of this compendium and its functions are built on the tidyverse language. The here package starts file pathing from the root folder of my scari package, allowing easier sharing. The rgbif package will be used to query records from GBIF, while the lydemapr package will be used to retrieve field survey data in the United States. GeoThinneR makes for easy spatial thinning of points and can employ a grid-based algorithm, so I will use this package.

Validation points- globally important viticultural regions

First, I will visualize the points we will use to validate our future MaxEnt models. These were retrieved by Huron et al, 2022 and represent 1,086 of the most important global viticultural regions. I will create a transformed version for use with the Behrmann CRS (ESRI:54017). Before doing that, I will deduplicate these records to ensure they are ready for use.

map_style <- list(
  xlab("UTM Easting"),
  ylab("UTM Northing"),
  theme_classic(),
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "lightblue2",
                                        colour = "lightblue2")
        ),
  scale_x_continuous(expand = c(0, 0)),
  scale_y_continuous(expand = c(0, 0)),
  viridis::scale_fill_viridis(option = "D"),
  coord_equal()
)

mypath <- file.path(here::here() %>% 
                     dirname(),
                   "maxent/historical_climate_rasters/chelsa2.1_30arcsec/v1_maxent_10km")

# background layer
global_bio2_df <- terra::rast(x = file.path(mypath, "bio2_1981-2010_global.asc")) %>%
  terra::as.data.frame(., xy = TRUE)

# wineries
load(file.path(here::here(), "data-raw", "wineries.rda"))

IVR_locations <- wineries %>%
  # tidy
  dplyr::select(-c(Latitude, Longitude)) %>%
  dplyr::filter(
    !is.na(x),
    !is.na(y)
    ) %>%
  dplyr::mutate(record_type = "IVR_location")

# remove duplicate coordinates
IVR_locations <- IVR_locations %>%
  CoordinateCleaner::clean_coordinates(
    x = .,
    lon = "x",
    lat = "y",
    species = "record_type",
    tests = "duplicates",
    value = "clean" # just return same df without duplicates
  ) %>%
  dplyr::select(., -record_type)

## Testing coordinate validity

## Flagged 0 records.

## Testing duplicates

## Flagged 7 records.

## Flagged 7 of 1086 records, EQ = 0.01.

The tidying removed 7 duplicate records, leaving 1,079 total records.

# I need to transform this dataset to the correct crs
IVR_locations_sv <- terra::vect(IVR_locations, geom = c("x", "y"), crs = "EPSG:4326") %>%
    terra::project(y = "ESRI:54017") %>% # convert to UTM 33N
# convert to geom, which gets coordinates of a spatVector
    terra::geom()
# convert back to data frame
IVR_locations_UTM <- terra::as.data.frame(IVR_locations_sv) %>%
    dplyr::select(-c(geom, part, hole)) %>%
    dplyr::rename("x_utm" = "x", "y_utm" = "y") 

# bind OG data
IVR_locations_UTM <- cbind(IVR_locations_UTM, IVR_locations) %>% 
  dplyr::select(-c(x, y)) %>% 
  dplyr::rename("x" = "x_utm", "y" = "y_utm") %>%
  # create ID col
  dplyr::mutate(ID = row_number()) %>%
  relocate(ID, x, y)

IVR_points_plot <- ggplot() +
    geom_raster(data = global_bio2_df, aes(x = x, y = y), fill = "azure1") +
    geom_point(data = IVR_locations_UTM, aes(x = x, y = y), fill = "purple3", shape = 21, size = 1, stroke = 0.3) +
    ggtitle("globally important viticultural regions") +
    map_style +
    theme(
      legend.position = "none",
      axis.title = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank()
      )

# save plot
ggsave(
  IVR_points_plot, 
  filename = file.path(here::here(), "vignette-outputs", "figures", "global_viticultural_regions.jpg"),
  height = 8, 
  width = 12,
  device = jpeg,
  dpi = "retina"
  )

# save transformation for future use
readr::write_rds(x = IVR_locations_UTM, file = file.path(here::here(), "data", "wineries_esri54017.rds"))

1. Retrieve data from GBIF and lydemapR and tidy

1.1- GBIF via rgbif

I will begin by retrieving the GBIF taxonomic ID for Lycorma delicatula.

slf_id <- rgbif::occ_search(scientificName = "Lycorma delicatula")[["data"]]

slf_id <- slf_id %>%
  dplyr::select(taxonKey) %>%
  dplyr::slice_head() %>%
  as.character()

I will also perform a few operations specific to rgbif. I will need to enter my user credentials for gbif; to keep these private, I will edit the .Renviron and call them from there. Once this chunk is run and the .Renviron pops up, enter your username, password and email credentials in the following format on lines 1-3:

GBIF_USER=” ”
GBIF_PWD=” ”
GBIF_EMAIL=” ”

Save and close the document.

# edit R environment for user credentials
usethis::edit_r_environ()

# get list of countries
rgbif::enumeration_country(curlopts = list())

Retrieve and save records

The retrieved ID is 5157899. This matches the ID for the SLF repository on GBIF.

Next, retrieve GBIF records using occ_download(). I will be editing the following GBIF constraints:

hasCoordinate: limit records to only those with coordinate data
year: we will retrieve records between 1981 and 2023. 1981 was chosen as the earliest date because it corresponds with the earliest climate data available.
occurrenceStatus: set to only presences, because absence data is not needed
basisOfRecord: here we include everything except fossil records, to ensure that records are within the target time period.
hasGeospatialIssue: exclude any records with an issue in the coordinate data
country: the country of the record

For the countries, we will reference the GBIF repo for SLF above, looking specifically at the available data tables and including countries hosting known populations. I list the counties inclusively, but I will intentionally leave out countries in which no known populations exist. For example, Mexico and Canada were excluded on the basis of expert opinion on the location of the suspected locations in GBIF. Taiwan was also excluded based on a paper by Lin and Liao- this paper performed a comprehensive survey for SLF in Taiwan to validate SLF presence records in online repositories and found that these records were likely false (Lin and Liao, 2024). Finally, I also also compared the map of points available from both GBIF and iNaturalist. GBIF contains most iNaturalist points, excluding a few protected under certain licenses which are not publicly available.

 # countries: include China, Japan, South Korea, North Korea, India, Vietnam, Bangladesh, USA
countries_iso <- c("CN", "JP", "KR", "KP", "IN", "VN", "BN", "US")

# initiate download
slf_gbif <- rgbif::occ_download(
  # general formatting
  type = "and",
  format = "SIMPLE_CSV",
  # inclusion rules
  pred("taxonKey", slf_id), # search by ID, not species name
  pred("hasCoordinate", TRUE),
  pred("hasGeospatialIssue", FALSE),
  pred("occurrenceStatus", "PRESENT"),
  pred_in("country", countries_iso),
  pred_gte("year", 1981), # records from 1981 onwards
  # exclusion rules
  pred_not(pred_in("basisOfRecord", "FOSSIL_SPECIMEN")) # exclude fossil records
)

# be sure to open slf_gbif and retrieve the download key

# check on download status using download key
rgbif::occ_download_wait('0000043-260623161305970')

The data have been pulled to the gbif server, but now we must download the data and edit it in a number of ways. I will use the download key for all of my functions from now on: key = '0000043-260623161305970'.

First, I will import the dataset to my working directory. I will also retrieve the data citation for this download. I will rename the raw download file. Finally, I will subset the desired columns and save them as a separate .csv for use in tidying. The raw query and the raw coordinates will be saved to the data-raw folder. I will read in the coordinate data saved there afterwards to prepare for data tidying.

# download the data as a zip to PC
slf_gbif_download <- rgbif::occ_download_get(
  key = '0000043-260623161305970',
  path = file.path(here::here(), "data-raw"),
  # overwrite = TRUE
  ) %>% 
  rgbif::occ_download_import()


# .zip file
# unzip
utils::unzip(
  zipfile = file.path(here::here(), "data-raw", "0000043-260623161305970.zip"),
  exdir =  file.path(here::here(), "data-raw")
  )
# rename .csv file
file.rename(
  from = file.path(here::here(), "data-raw", "0000043-260623161305970.csv"),
  to = file.path(here::here(), "data-raw", paste0("slf_gbif_raw_query_", format(Sys.Date(), "%Y-%m-%d"), ".csv"))
  )
# delete original zip
file.remove(file.path(here::here(), "data-raw", "0000043-260623161305970.zip"))



# simplify download into more useable .csv
slf_gbif_download_simple <- slf_gbif_download %>%
  mutate(prov = "gbif") %>%
  dplyr::select(c("species", "decimalLongitude", "decimalLatitude", "countryCode", "stateProvince", "prov", "year", "gbifID")) %>% # select desired columns
  # rename for consistency with rest of vignette
  rename("name" = "species",
         "longitude" = "decimalLongitude", 
         "latitude" = "decimalLatitude",
         "key" = "gbifID"
  )

# save .csv both outputs for use
write_csv(x = slf_gbif_download_simple, 
          file = file.path(here::here(), "data-raw", paste0("slf_gbif_raw_coords_", format(Sys.Date(), "%Y-%m-%d"), ".csv"))
)

rgbif::gbif_citation('0000043-260623161305970')

Here is the data citation: “GBIF Occurrence Download https://doi.org/10.15468/dl.wfpdrz Accessed from R via rgbif (https://github.com/ropensci/rgbif) on 2026-06-23”

There are 59,999 total records from GBIF.

Coordinate veracity

First, we want to check the coordinate veracity using the CoordinateCleaner package. We will be more picky and apply more tests to the GBIF data points than those from the other datasets because they are usually less accurate than field survey data.

valid and likely coordinates
points set to the centroid of a country or province
coordinates that fall into the ocean
duplicates
points with equal lat / lon
points near the GBIF HQ
points near biodiv institutions
points near whole lat / lon numbers

# read in data
slf_gbif_coords1 <- readr::read_csv(file = file.path(here::here(), "data-raw", "slf_gbif_raw_coords_2026-06-23.csv"))

# set random seed
set.seed(999)

# coordinate cleaner
slf_gbif_report <- CoordinateCleaner::clean_coordinates(
  x = slf_gbif_coords1, 
  lon = "longitude",
  lat = "latitude",
  species = "name",
  tests = c("centroids", "duplicates", "equal", "gbif", "institutions", "seas", "zeros"),
  # test details
  seas_scale = 110,
  centroids_detail = "both",
  verbose = TRUE
) %>%
  as.data.frame()

It seems like the most records were removed for duplicate records, followed by for records flagged as being in the ocean or other bodies of water. We also did find a few records that were attributed to biodiversity institutions and country / province centroids.

Next, I will filter out problematic points and ensure the species naming is consistent. We are removing 6261 records in total.

# remove flagged records
slf_gbif_coords1 <- slf_gbif_report %>%
  dplyr::filter(.summary == TRUE) # FALSE records are potentially problematic

# check species name
unique(slf_gbif_coords1$name)

# there is only one naming convention

Next, we will use GeoThinneR::thin_points to spatially thin the occurrence data. The GeoThinneR package was chosen to perform the initial thinning because it preserves non-coordinate data in its output. So, we can save a spatially thinned version of the data with other information (such as unique keys or dates) that might be useful later.

In GeoThinneR::thin_points(), The data will be rarefied to 10km over 10 passes. This will be done once more on the final dataset. Spatial thinning should reduce autocorrelation and sampling bias. I have programmed this function to rarefy at 10km, which is the current resolution of the climate data we used from CHELSA.GeoThinneR has a grid function that will thin species data based on a predetermined raster, so we will use the CHELSA rasters we tidied and aggregated in the last vignette. These rasters originally came at a resolution of 30 arc-seconds (roughly 1km at the equator) but we aggregated them at 10km by taking the mean value per cell. We will also perform this process separately for different countries. I am using a special reference layer I created that is identical to the rasters that will be used for the MaxEnt models, except that it is projected into the EPSG:4326 coordinate reference system (CRS), instead of the ESRI:54107 crs.

geothinner_raster <- terra::rast(x = file.path(mypath, "geothinning_ref_layer_epsg4326_bio2.asc"))

slf_gbif_coords2 <- GeoThinneR::thin_points(
  data = slf_gbif_coords1,
  long_col = "longitude",
  lat_col = "latitude",
  method = "grid", # use a grid based thinning method
  raster_obj = geothinner_raster, # the aggregated raster to use
  group_col = "countryCode", # additionally thin separately based on country
  trials = 10, # number of passes
  seed = 997, # set reproducible seed
  verbose = TRUE
  ) 

View(slf_gbif_coords2[[1]])

# get rid of unnecessary columns used for thinning
slf_gbif_coords2 <- slf_gbif_coords2[[1]] %>% dplyr::select(name:key)

Spatial thinning left 3020 of 53,738 total records. The last action we will perform is to filter the data by state. We have compiled a list of US States with verified, established populations of Lycorma Delicatula. This was provided and validated by the USDA.

# object containing state names
USDA_states <- c("Connecticut", "District of Columbia", "Delaware", "Georgia", "Illinois", "Indiana", "Kentucky", "Maryland", "Massachusetts", "Michigan", "New Jersey", "New York", "North Carolina", "Ohio", "Pennsylvania", "Rhode Island", "South Carolina", "Tennessee", "Virginia", "West Virginia")

# next, check for naming duplicates
slf_gbif_coords2 %>% 
  dplyr::filter(.$countryCode == "US") %>%
  dplyr::select("stateProvince") %>%
  as.data.frame() %>%
  unique()

There are 29 unique states where SLF is found within the GBIF database (including DC). However, we will only keep records from the 19 states with verified populations, including DC.

slf_gbif_coords2 <- slf_gbif_coords2 %>%
  dplyr::filter(.$countryCode != "US" | .$stateProvince %in% USDA_states)  # keep only records which are not in the USA, or records which are in the USA and are within the USDA list of states. 

# run this code again to ensure filtering worked
slf_gbif_coords2 %>% 
  dplyr::filter(.$countryCode == "US") %>%
  dplyr::select("stateProvince") %>%
  as.data.frame() %>%
  unique()

# it worked!

We will save the results of the cleaning we just performed so they can be referenced or called later in our analysis. Later in this vignette, we will be combining these records with those from other data sources into a final dataset that is prepared for MaxEnt.

if(TRUE){
  
  write_csv(slf_gbif_coords2, file = file.path(here::here(), "vignette-outputs", "data-tables", "slf_gbif_cleaned_coords_2026-06-23.csv"))
  
  # also save as a .rds file
  write_rds(slf_gbif_coords2, file = file.path(here::here(), "data", "slf_gbif_cleaned_coords_2026-06-23.rds"))
  
}

1.2- LydemapR

Retrieve and save lydemapR records

Next, we will retrieve data from the lydemapR package and save it as a .csv. I will repeat the spatial thinning and cleaning process outlined above.

slf_lyde <- lydemapr::lyde

write_csv(x = slf_lyde, 
          file = file.path(here::here(), "data-raw", paste0("slf_lyde_raw_coords_", format(Sys.Date(), "%Y-%m-%d"), ".csv")))

The next chunk was an import from an unpublished version of the lydemapr dataset. This chunk is not being used as of scari 1.0.0

if (FALSE) {
  
  slf_lyde_raw <- readr::read_csv(file = file.path(here::here(), "data-raw", ""))
  
  # overwrite save original version
  write_csv(x = slf_lyde_raw, file = file.path(here::here(), "data-raw", ""))
  
}

The raw dataset is composed of over 1 million records of SLF in the USA alone. Most of these records are concentrated in the mid-atlantic region where the invasion front is progressing. This dataset will also need a run of spatial thinning. First, we will need to narrow these records to fit our needs. At the time of writing this, I am only interested in records that were collected during field surveys, but this can be tuned to the needs of the analysis. So, I will set the collection method equal to the value of “field_survey/management”. Obviously, I am also interested in presence records only.

Lastly, I will pull records only from areas where SLF have established a population. LydemapR defines an established population as either 2+ adults or the presence of 1 egg mass.

# read in data
slf_lyde <- readr::read_csv(file = file.path(here::here(), "data-raw", "slf_lyde_raw_coords_2026-07-30.csv"))

# unique values for collection method
unique(slf_lyde$collection_method)

slf_lyde1 <- slf_lyde %>%
  dplyr::filter(lyde_present == "TRUE", # only select presences
         collection_method == "field_survey/management", # collection method = "field_survey/management"
         lyde_established == "TRUE") %>% # only select areas records that are from established populations
  # add species column
  dplyr::mutate(species = "Lycorma delicatula")

The data that fit our needs are about 81,700 records total.

Coordinate veracity

Here we will run through the same data cleaning process that we performed for the GBIF data.

set.seed(996)

# coordinate cleaner
slf_lyde_report <- CoordinateCleaner::clean_coordinates(
  x = slf_lyde1, 
  lon = "longitude",
  lat = "latitude",
  species = "species",
  tests = c("centroids", "duplicates", "equal", "gbif", "institutions", "seas", "zeros"),
  # test details
  seas_scale = 110,
  centroids_detail = "both"
) %>%
  as.data.frame()

# remove flagged records
slf_lyde1 <- slf_lyde_report %>%
  dplyr::filter(.summary == TRUE) # FALSE records are potentially problematic

Again, the majority were duplicates, followed by those that fell into bodies of water. This removal left about 12,700 records.

slf_lyde2 <- GeoThinneR::thin_points(
  data = slf_lyde1,
  long_col = "longitude",
  lat_col = "latitude",
  method = "grid", # use a grid based thinning method
  raster_obj = geothinner_raster, # the aggregated raster to use
  trials = 10, # number of passes
  seed = 995, # set reproducible seed
  verbose = TRUE
  ) 

# get rid of unnecessary columns used for thinning
slf_lyde2 <- slf_lyde2[[1]] %>% dplyr::select(source:species)

Spatial thinning for the lyde dataset left 1,673 final points that fit our needs. Previous usage of this package left 1371 points out of 39,000 after the initial selection of data (only established populations, etc.).

write_csv(slf_lyde2, file = file.path(here::here(), "vignette-outputs", "data-tables", "slf_lyde_cleaned_coords_2026-07-30.csv"))

# also save as a .rds file
write_rds(slf_lyde2, file = file.path(here::here(), "data", "slf_lyde_cleaned_coords_2026-07-30.rds"))

1.3- Visualize SLF distributions

Visualize results of spatial thinning

First, I will ensure that the spatial thinning was correct and efficient. We want to visualize the difference between the raw coordinates that were pulled from GBIF and the cleaned, spatially thinned version to ensure proper geospatial coverage and the success of the spatial thinning runs.

# raw data
gbif_raw <- readr::read_csv(file = file.path(here::here(), "data-raw", "slf_gbif_raw_coords_2026-06-23.csv"))

# thinned data
gbif_thinned <- readr::read_csv(file = file.path(here::here(), "vignette-outputs", "data-tables", "slf_gbif_cleaned_coords_2026-06-23.csv"))

# plot world map
map_gbif_thinned_NAmerica <- ggplot() +
  # basemap
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = gbif_raw, aes(x = longitude, y = latitude), color = "blue", size = 2) +
  geom_point(data = gbif_thinned, aes(x = longitude, y = latitude), color = "darkorange", shape = 2) +
  #coord_sf(xlim = c(-164.5, 163.5), ylim = c(-55, 85)) +
  coord_sf(xlim = c(-133.593750, -52.294922), ylim = c(25.085599, 55.304138)) +
  ggtitle("SLF GBIF records in N America- \n raw (blue) vs thinned (darkorange)") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw() +
  theme(axis.title = element_blank())

# plot map of Asia
map_gbif_thinned_Asia <- ggplot() +
  # basemao
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = gbif_raw, aes(x = longitude, y = latitude), color = "blue", size = 2) +
  geom_point(data = gbif_thinned, aes(x = longitude, y = latitude), color = "darkorange", shape = 2) +
  coord_sf(xlim = c(100, 140), ylim = c(20, 50)) +
  ggtitle("SLF GBIF records in Asia-\nraw (blue) vs thinned (darkorange)") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw() +
  theme(axis.title = element_blank())
  
# patchwork display of plots
map_gbif_thinned_NAmerica + map_gbif_thinned_Asia +
  plot_layout(nrow = 2)

From the mapping, it seems that the spatial thinning reduced the GBIF records, while preserving the same spatial extent. We also zoom in on the invaded range, and indeed every blue point is paired with at least 1 dark orange point. One exception is a blue point located near Jeju island, south of the mainland, but this is likely due to coastline inconsistencies (normal when aligning data points to a raster). So, our goal was met!

Visualize results of excluding records not representing established populations

For the Lydemapr package data, I chose to include only established populations so that the data would be more reliable and less repetitive (the sightings not representing established populations might represent regulatory incidents, and definitely represent the GBIF data already retrieved because these are integrated into LydemapR). However, I want to ensure there is no unreasonable loss of geographic area based on this, so we will map the points.

# raw data
lyde_raw <- readr::read_csv(file = file.path(here::here(), "data-raw", "slf_lyde_raw_coords_2026-07-30.csv"))

# ensure data points represent presences
lyde_presences <- lyde_raw %>%
  filter(lyde_present == "TRUE") 

# include only established records 
lyde_establishments <- lyde_raw %>%
  filter(
    lyde_present == "TRUE",
    lyde_established == "TRUE"
    ) %>% 
  mutate(species = "Lycorma delicatula")

# plot map of NAmerica
map_lyde_NAmerica <- ggplot() +
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = lyde_presences, aes(x = longitude, y = latitude), color = "blue", size = 2) +
  geom_point(data = lyde_establishments, aes(x = longitude, y = latitude), color = "darkorange", shape = 2) +
  coord_sf(xlim = c(-133.593750, -52.294922), ylim = c(25.085599, 55.304138)) +
  ggtitle("SLF lydemapR records- \nall records (blue)\nvs established populations (dark orange) for N America") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw()
  
# plot map including only established records
map_lyde_NEast <- ggplot() +
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = lyde_presences, aes(x = longitude, y = latitude), color = "blue", size = 2) +
  geom_point(data = lyde_establishments, aes(x = longitude, y = latitude), color = "darkorange", shape = 2) +
  coord_sf(xlim = c(-90, -70), ylim = c(35, 45)) +
  ggtitle("SLF lydemapR records- \nall records (blue)\nvs established populations (dark orange) for the eastern USA") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw()

# patchwork display of plots
 map_lyde_NAmerica + map_lyde_NEast +
  plot_layout(nrow = 2)

Indeed, we did lose some geographical area in the USA, but these are likely not established populations or are areas where SLF could not reasonably establish a population (expert opinion).

2. Combine records

2.1- SLF data from published sources

In this step, I will combine the records obtained from GBIF and LydemapR with data taken from the literature. I will begin by loading in records taken from the published literature. There are a total of 217 records from China, southeast Asia, Japan, SK and the USA. Most of these records were obtained from genetic studies and represent samples taken from established populations of SLF. Data is scarce for established populations of SLF within its native range (China and southeast Asia), so these records are especially important. These records were taken from both peer-reviewed literature and natural history notes.

The key column was created to house the unique key ID from each paper’s samples. If a paper did not specify one, I created a unique key using a the following convention: Author last name (1st 3 letters), year of publication, and a number to represent the sample.

I had to transform many records from other latitude/longitude conventions using a coordinate converter. I converted to the decimal degree (DD) format. I used a few sites: - epsg.io: https://epsg.io/transform#s_srs=4326&t_srs=54017&x=NaN&y=NaN - UMinn PCG: https://applications.pgc.umn.edu/convert/

# read in .csv of papers
slf_published_papers <- readr::read_csv(file.path(here::here(), "data-raw", "slf_publishedOccurrenceRecords_papers.csv"))

# kable table

slf_published_papers %>%
  dplyr::select(!notes) %>%
  knitr::kable(format = "pipe")

author	year	record_spatial_extent	DOI	key_ID
Fennah	1953	CN	NA	Fen_1953_##
Han et al	2008	SK	10.1111/j.1748-5967.2008.00188.x	Han_2008_##
Kim et al	2013	CN, JP, SK	10.1016/j.aspen.2013.07.003	included
Yang et al	2015	CN	10.11865/zs.20150305	Yang_2015_##
Zhang et al	2019	CN	10.3390/insects10100312	included
Xin et al	2020	CN	10.1093/ee/nvaa137	Xin_2020_##
Du et al	2021	CN, JP, SK, US	10.1111/eva.13170	included
Kim et al	2021	CN, JP, SK	10.3390/insects12060539	included
Manzoor et al	2021	CN	10.1111/aab.12674	Man_2021_##
Nakashita et al	2022	CN, JP, SK, US	10.1038/s41598-022-05541-z	included
Suzuki	2023	JP	NA	Suz_2023_##
Kamiyama et al	2024	JP	10.1002/1438-390X.12203	included
Liu et al	2024	US	10.1016/j.cris.2024.100078	Liu_2024_##
Kamiyama et al	2026	CN, JP, SK, US	10.1002/ecs2.70630	Kam_2026_##

slf_published_papers

## # A tibble: 14 × 6
##    author           year record_spatial_extent DOI                  key_ID notes
##    <chr>           <dbl> <chr>                 <chr>                <chr>  <chr>
##  1 Fennah           1953 CN                    NA                   Fen_1… repr…
##  2 Han et al        2008 SK                    10.1111/j.1748-5967… Han_2… NA   
##  3 Kim et al        2013 CN, JP, SK            10.1016/j.aspen.201… inclu… NA   
##  4 Yang et al       2015 CN                    10.11865/zs.20150305 Yang_… repr…
##  5 Zhang et al      2019 CN                    10.3390/insects1010… inclu… NA   
##  6 Xin et al        2020 CN                    10.1093/ee/nvaa137   Xin_2… NA   
##  7 Du et al         2021 CN, JP, SK, US        10.1111/eva.13170    inclu… NA   
##  8 Kim et al        2021 CN, JP, SK            10.3390/insects1206… inclu… NA   
##  9 Manzoor et al    2021 CN                    10.1111/aab.12674    Man_2… repr…
## 10 Nakashita et al  2022 CN, JP, SK, US        10.1038/s41598-022-… inclu… NA   
## 11 Suzuki           2023 JP                    NA                   Suz_2… NA   
## 12 Kamiyama et al   2024 JP                    10.1002/1438-390X.1… inclu… NA   
## 13 Liu et al        2024 US                    10.1016/j.cris.2024… Liu_2… coor…
## 14 Kamiyama et al   2026 CN, JP, SK, US        10.1002/ecs2.70630   Kam_2… coor…

slf_published <- readr::read_csv(file.path(here::here(), "data-raw", "slf_publishedOccurrenceRecords_2026-06-25.csv"))

## Rows: 217 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): name, key, country, stateProvince, publishingArticle, accessionNum,...
## dbl (3): longitude, latitude, year
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# save as .rds before proceeding
write_rds(x = slf_published, file = file.path(here::here(), "data", "slf_publishedOccurrenceRecords_2026-06-25.rds"))

map_published_Asia <- ggplot() +
  # basemap
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = slf_published, aes(x = longitude, y = latitude), color = "darkorange") +
  coord_sf(xlim = c(100, 140), ylim = c(10, 45)) +
  ggtitle("SLF data from published sources \n in China") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw() +
  theme(axis.title = element_blank())

# map alongside previous map of GBIF points

# bounding box coords found at: 
# http://bboxfinder.com/#0.000000,0.000000,0.000000,0.000000

# patchwork
map_published_Asia + map_gbif_thinned_Asia +
  plot_layout(nrow = 2)

The map above shows that the occurrence data from the literature provide some coverage of the native range that we might not otherwise get especially in the west and south of china.

map_published_NAmerica <- ggplot() +
  # basemap
  geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
  geom_point(data = slf_published, aes(x = longitude, y = latitude), color = "darkorange") +
  coord_sf(xlim = c(-133.593750, -52.294922), ylim = c(25.085599, 55.304138)) +
  ggtitle("SLF data from published sources \n in N America") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  theme_bw() +
  theme(axis.title = element_blank())

map_published_NAmerica

We also see that there were very few records in the USA- these likely did not change our geographic extent much.

2.2- Tidy and join datasets

Now, I will load in the cleaned coordinates from GBIF and LydemapR that we produced earlier.

slf_gbif_coords3 <- readr::read_csv(file.path(here::here(), "vignette-outputs", "data-tables", "slf_gbif_cleaned_coords_2026-06-23.csv"))

slf_lyde3 <- readr::read_csv(file.path(here::here(), "vignette-outputs", "data-tables", "slf_lyde_cleaned_coords_2026-07-30.csv"))

First, I will tidy the data for joining. We will keep only the species names, coordinates and unique key ID. We will also add a column that states which data source each point is from.

# tidy gbif data
slf_gbif_coords3 %<>%
  dplyr::select(name:latitude, key, prov) %>%
  rename("data_source" = "prov",
         "species" = "name")

# tidy lyde data
slf_lyde3 %<>%
  dplyr::select(species, longitude, latitude, pointID) %>%
  rename("key" = "pointID") %>%
  mutate(data_source = "lyde")
  
# published data
slf_published %<>%
  dplyr::select(name:key, publishingArticle) %>%
  rename("species" = "name",
         "data_source" = "publishingArticle")

# the publishing data source column needs to be tidied further. I will take out commas and substitute spaces for underscores
slf_published$data_source <- gsub(pattern = " ", replacement = "_", x = slf_published$data_source) 
slf_published$data_source <- gsub(pattern = ",", replacement = "", x = slf_published$data_source)

# Finally, use head() to check coltypes are the same
head(slf_lyde3)
head(slf_gbif_coords3)
head(slf_published)

# we see that the key column in the GBIF dataset is a double, while the key columns in the other 2 are characters. We will change this:
slf_gbif_coords3$key <-  as.character(slf_gbif_coords3$key)

head(slf_gbif_coords3)

# the conversion worked

3. Final Data Tidying

Now the data are in a format that is easier to join. We will join the data and save the cleaned coordinates.

slf_all_coords <- slf_gbif_coords3 %>%
  full_join(., slf_lyde3) %>%
  full_join(., slf_published)

write_csv(x = slf_all_coords, file = file.path(here::here(), "vignette-outputs", "data-tables", "slf_all_coords_2026-07-30.csv"))

3.1- Final spatial thinning

The last step of tidying these data is to perform a second round of spatial thinning. We need to do this again because after joining 3 different datasets together, some of the points may be within 10km of each other. The data will be rarefied to 10km over 10 passes. This will be done twice more on the final dataset. Spatial thinning should reduce autocorrelation and sampling bias. Again, we will thin to a minimum of 10km distance between points. We will repeat the process above of using CoordinateCleaner and GeoThinneR.

Both files will be written to the vignettes-outputs folder, which holds intermediate data objects that are not the final versions of the data.

slf_all_coords1 <- GeoThinneR::thin_points(
  data = slf_all_coords,
  long_col = "longitude",
  lat_col = "latitude",
  method = "grid", # use a grid based thinning method
  raster_obj = geothinner_raster, # the aggregated raster to use
  trials = 10, # number of passes
  seed = 995, # set reproducible seed
  verbose = TRUE
  ) 

slf_all_coords1 <- slf_all_coords1[[1]]

This thinning method left 3422 of 4892 records.

# coordinate cleaner
slf_all_coords_report <- CoordinateCleaner::clean_coordinates(
  x = slf_all_coords1, 
  lon = "longitude",
  lat = "latitude",
  species = "species",
  tests = c("centroids", "duplicates", "equal", "gbif", "institutions", "seas", "zeros")
) %>%
  as.data.frame()

# remove flagged records
slf_all_coords2 <- slf_all_coords_report %>%
  dplyr::filter(.summary == TRUE) %>% # FALSE records are potentially problematic
  dplyr::select(species:latitude)

This thinning method removed an additional ~80 records, reducing the final number to 3345.

Finally, I will re-run the GeoThinneR method, additionally using a raw distance metric.

slf_all_coords3 <- GeoThinneR::thin_points(
  data = slf_all_coords2,
  long_col = "longitude",
  lat_col = "latitude",
  method = "brute_force",
  euclidean = FALSE, # use the haversine distance metric because these coordinates are projected to a lat/long system
  thin_dist = 10,
  seed = 995, # set reproducible seed
  verbose = TRUE
  ) 

slf_all_coords3 <- slf_all_coords3[[1]]

This method removed a lot of coordinates- reducing the final count from 3345 to 1967 records.

3.2- Save output for MaxEnt

We will finally convert the slf coordinates to the format that will be used by MaxEnt and save it. That is, the first column “species” must be the species name, the second column “x” must be the longitude, and the third column “y” must be the latitude. We will also need to transform the coordinates to the ESRI:54017 coordinate reference system (CRS) we used in the previous vignette.

# save EPSG4326 version before transformation
write_csv(x = slf_all_coords3, file = file.path(here::here(), "vignette-outputs", "data-tables", paste0("slf_all_coords_final_2026-07-30_epsg4326.csv")))

# rename columns
slf_all_coords4 <- slf_all_coords3 %>%
  rename("x" = "longitude",
         "y" = "latitude")

# transform CRS of downloaded data
slf_all_coords4 <- terra::vect(slf_all_coords4, geom = c("x", "y"), crs = "EPSG:4326") %>%
    terra::project(y = "ESRI:54017") %>% # convert to UTM 33N
# convert to geom, which gets coordinates of a spatVector
    terra::geom()
# convert back to data frame
slf_all_coords5 <- terra::as.data.frame(slf_all_coords4) %>%
  dplyr::select(-c(geom, part, hole)) %>%
  # join original dataset back for data
  cbind(., slf_all_coords3) %>%
  # remove lat long and other edits
  dplyr::select(-c(longitude, latitude)) %>%
  dplyr::relocate(species)

# save
write_csv(x = slf_all_coords5, file = file.path(here::here(), "vignette-outputs", "data-tables", paste0("slf_all_coords_final_2026-07-30.csv")))

# also save as a .rds file
write_rds(slf_all_coords5, file = file.path(here::here(), "data", paste0("slf_all_coords_final_2026-07-30.rds")))

slf_all_coords_df <- readr::read_csv(file = file.path(here::here(), "vignette-outputs", "data-tables", "slf_all_coords_final_2026-07-30_epsg4326.csv"))

## Rows: 1967 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): species
## dbl (2): longitude, latitude
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

 map_all <- ggplot() +
    # basemap
    geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
    geom_point(data = slf_all_coords_df, aes(x = longitude, y = latitude), color = "darkorange", size = 1) +
    coord_sf(xlim = c(-164.5, 163.5), ylim = c(-55, 85)) +
    ggtitle("All SLF presences") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank()) +
    theme_bw() +
    theme(axis.title = element_blank())

 map_all_NAmerica <- ggplot() +
    # basemap
    geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
    geom_point(data = slf_all_coords_df, aes(x = longitude, y = latitude), color = "darkorange", size = 1) +
    coord_sf(xlim = c(-133.593750, -52.294922), ylim = c(25.085599, 55.304138)) +
    ggtitle("SLF presences N America") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank()) +
    theme_bw() +
    theme(axis.title = element_blank())
 
  map_all_Asia <- ggplot() +
    # basemap
    geom_sf(data = countries_sf, fill = NA, color = "black", lwd = 0.15) +
    geom_point(data = slf_all_coords_df, aes(x = longitude, y = latitude), color = "darkorange", size = 1) +
    coord_sf(xlim = c(68.906250, 152.534180), ylim = c(8.928487, 45.920587)) +
    ggtitle("SLF presences Asia") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank()) +
    theme_bw() +
    theme(axis.title = element_blank())

# patchwork
map_all_NAmerica + map_all_Asia +
plot_layout(nrow = 2)

map_all

From the map, we can see the final points that were selected. We have a good spread across the native range and better spread than before in the invaded range. This version of the total dataset included newer data from lydemapr, so we now have better coverage in the invaded range. I will also check the number of points in the native vs the invaded range.

SLF_by_longitude <- ggplot(data = slf_all_coords_df) +
  geom_histogram(aes(x = longitude)) +
  ggtitle("count of SLF presence records invaded vs native range") +
  labs(caption = "the longitudinal range between the two dashed lines represents the native range for SLF") +
  xlab("longitude") +
  ylab("count of SLF presences") +
  geom_vline(xintercept = 73.37, linetype = "dashed") +
  geom_vline(xintercept = 124.06, linetype = "dashed") +
  theme_minimal()

SLF_by_longitude

# save output
ggsave(
  SLF_by_longitude, filename = file.path(here::here(), "vignette-outputs", "figures", "SLF_all_coords_final_by_longitude.jpg"),
  height = 8, 
  width = 10,
  device = jpeg,
  dpi = 500
  )

The native range includes longitudes between 35 and 105. We can see that most of the SLF presence data are from the invaded ranges of Japan, South Korea or the United States. This is a weakness in the data available and we have given our best effort to account for it by including presence data for China from the literature.

We now have half of the data we need to perform SDM! In the next vignette, we will download historical and projected climate and human impact data, which are the basis of our methods for predicting suitability for Lycorma delicatula.

References

De Bona, S., Barringer, L., Kurtz, P., Losiewicz, J., Parra, G. R., & Helmus, M. R. (2023). lydemapr: An R package to track the spread of the invasive spotted lanternfly (Lycorma delicatula, White 1845) (Hemiptera, Fulgoridae) in the United States. NeoBiota, 86, 151–168. https://doi.org/10.3897/neobiota.86.101471
GBIF Occurrence Download https://doi.org/10.15468/dl.59gndc Accessed from R via rgbif (https://github.com/ropensci/rgbif) on 2024-08-05
Huron, N. A., Behm, J. E., & Helmus, M. R. (2022). Paninvasion severity assessment of a U.S. grape pest to disrupt the global wine market. Communications Biology, 5(1), 655. https://doi.org/10.1038/s42003-022-03580-w
Lin, Y.-S., & Liao, J.-R. (2024). Multifaceted Investigation into the Absence and Potential Invasion of Spotted Lanternfly (Lycorma delicatula) in Taiwan. Research Square. https://doi.org/10.21203/rs.3.rs-4832573/v1
Mestre-Tomás, J. (2024). GeoThinneR: An R package for simple spatial thinning methods in ecological and spatial analysis. R package version 1.0.0, https://github.com/jmestret/GeoThinneR
epsg.io, powered by MapTiler. Transform coordinates - GPS online converter [Internet]. Available from: https://epsg.io/transform#s_srs=4326&t_srs=54017&x=NaN&y=NaN
Coordinate Converter – Polar Geospatial Center [Internet]. Available from: https://applications.pgc.umn.edu/convert/

Samuel M. Owens1

2026-07-30

Overview

Setup

Validation points- globally important viticultural regions

1. Retrieve data from GBIF and lydemapR and tidy

1.1- GBIF via rgbif

Retrieve and save records

Coordinate veracity

1.2- LydemapR

Retrieve and save lydemapR records

Coordinate veracity

1.3- Visualize SLF distributions

Visualize results of spatial thinning

Visualize results of excluding records not representing established populations

2. Combine records

2.1- SLF data from published sources

2.2- Tidy and join datasets

3. Final Data Tidying

3.1- Final spatial thinning

3.2- Save output for MaxEnt

References

Samuel M. Owens¹