Pennsylvania AFS Workshop

The Biology

Aaron Henning’s PA AFS Presentation from the Meeting on 2021-02-11 gives excellent background on the background about the map we’re about to make. Click here to download a pdf of Aaron's presentation

The Data Science

Download the R script we will use to make our maps first thing. Click here to download the script

Using R as GIS

Tabular Data

Often we just want to project point data from a simple tabular format - like a spreadsheet or rds file. This is super simple, we just import the tabular data into R and then use the x and y coordinates (latitude and longitude stored as decimal degrees) to project the points.

Geospatial Data

There are a number of packages, such as sf that we will use today, that can handle geospatial data such as GIS shapefiles. Geoprocessing operations, such as intersects, clips, buffers, etc. are accomplished with coding, rather than traditional point and click geoprocessing tools in ESRI applications.

We will first read in various types of spatial data using the sf package to deomonstrate this utility, and eventually build our map.

If you have a shapefile saved in your directory - or download one from github into your directory as we will do, use st_read() to import it into R. Then check the projection with st_crs(), and transform it to WGS 84 using st_transform() since that is the projection we need for our leaflet map.

nhd <- st_read('Flowlines.shp', stringsAsFactors = F) %>% # read in shapefile
  select(COMID, Resolution, GNIS_ID, GNIS_NAME, LENGTHKM, REACHCODE, FTYPE, StreamOrde, geometry) # select only fields we want

## Reading layer `Flowlines' from data source `C:\R\trainings\2021-02-12_PAAFS_leaflet\fromGH\Rmd\Flowlines.shp' using driver `ESRI Shapefile'
## Simple feature collection with 303 features and 133 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -76.91974 ymin: 39.54369 xmax: -75.95318 ymax: 40.30795
## CRS:            4269

## select: dropped 125 variables (FDATE, WBAREACOMI, FCODE, StreamLeve, StreamCalc, …)

st_crs(nhd) # extract coordinate system (CS) information from an sf

## Coordinate Reference System:
##   User input: 4269 
##   wkt:
## GEOGCS["GCS_North_American_1983",
##     DATUM["North_American_Datum_1983",
##         SPHEROID["GRS_1980",6378137.0,298.257222101]],
##     PRIMEM["Greenwich",0.0],
##     UNIT["Degree",0.0174532925199433],
##     AUTHORITY["EPSG","4269"]]

nhd <- st_transform(nhd, "+proj=longlat +datum=WGS84") # transform projection to be WGS84 (same as other leaflet layers)

This routine will work for point, line, and polygon shapefiles!

Spatial Data Packages

There are also packages that make spatial data available, such as USAboundaries, which contains state, county, and municipal boundaries across the entire US. These packages are convenient since they don’t require us to import a shapefile, but they are sometimes prone to version issues that can cause error messages. For our example, we will import the counties in PA and MD in our study area.

Here we have saved these county boundaries to an object called edna_co.

edna_co <- USAboundaries::us_counties(resolution = "high", states = c("PA", "MD")) %>% # PA & MD
  filter(name %in% c('York', 'Lancaster', 'Harford', 'Cecil')) # specify 4 counties

Now we can plot the counties to visualize

plot(edna_co$geometry)

The Pipe! `%>%`

Pipes are a powerful tool for clearly expressing a sequence of multiple operations. The pipe basically means ‘then’. It allows many operations to be linked together and prevents the need for multiple, intermediate dataframes on the way to your final product. You can easily include a pipe in your code by holding Ctrl + Shift + M. You can also view many additional shortcuts by holding Alt + Shift + K.

Additional DS Skills

In addition to learning about leaflet maps, you’re going to hopefully pick up some additional skills in R today.

Sourcing external data into your R environment
- Since I can’t physically walk around and hand everyone a flash drive, I’ve pushed all of the data files we need today to github. This has additional benefit when it comes to version control, collaboration with colleagues, and web publishing.
- Check out the repository for this workshop to see where all our files are sourced.
Data manipulation skills
- We’re using tidyverse functions to get our geospatial data into shape. Also using the %>% to link our code together and write clean code
Creating interactive, shareable documents.
- Today we will produce a fully interactive html file that we can publish as a website and share widely.

Introduction to leaflet

Leaflet is a popular open-source library used by news and social media organizations, GIS specialists, and scientists. There is a ton of helpful material online that can assist basic to advanced map creation. Leaflet can accomodate any of the point, line, or polygon spatial data described above, making it extremely flexible.

Leaflet is powerful because you can display spatial data in an interactive format, which is espcially useful for complex datasets like the NSH-BCAT project that has multiple species and sampling events. Multiple static maps would be required, and even then would provide a fraction of the information that a leaflet map is able to provide.

Here is an extremley helpful website that is close to a one-stop-shop for leaflet

In the simplest sense, we can add a dot to a map and make it interactive. This map is interactive, but simple. Once you’ve installed and loaded leaflet, you can run this code to produce a map of the location of the first snakehead captured for this project.

leaflet() %>% # Create a map widget by calling leaflet()
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers( # Add layers (i.e., features) to the map by using layer functions 
    lng=-76.20638, lat=39.61467 , # specify lat/long
    popup='<a href = "https://nas.er.usgs.gov/queries/SpecimenViewer.aspx?SpecimenID=1633025"> USGS NAS Report </a>') # add info in a popup box - this is a hyperlink to the USGS NAS report of the capture

This is a start, but pretty useless. Let’s get into the weeds to learn how to crunch real spatial data to make a polished map filled with useful information.

After initializing the map with leaflet() we can start by adding multiple layers to our map.

First add basemaps (tiles) to your map with addTiles()
- There are many available via the leaflet.extras package. Here is a great viewer to help you pick.
- You can also add multiple tiles to each map you create.
I like to add some fancy elements from leaflet.extras before getting into the geospatial layers
- addResetMapButton() provides a reset button to reset the map to the max extent
- addSearchOSM() provides a search box to allow users to search for towns, coordinates, etc.
Now it’s time to get into the meat and potatoes of the map and add geospatial layers
- Add point, line, and/or polygon data with addCircleMarkers(), addPolylines(), and addPolygons(), respectively
- You can also get creative with other marker/icon options
- The terms vary by data type here, but we will go into detail in the code
- We will also make informative pop-up box appear on a click of the point
- The group term will allow layers to be turned on/off and overlaid, so pay attention here!
Once our geospatial layers are added, we’ll add legends using addLegend()
- We have a lot of options, from categorical bins using colorFactor() to continuous color ramps using colorNumeric()
- We have to create our legend prior to initiating the map!
- Pro-tip: include a group term in your addLegend() function and the legend will turn on/off with the layer!
Add a scale bar using addScaleBar() because why not?
Use addLayersControl() to manage multiple basemaps, geospatial layers, etc.
The map will initialize at the maximum extent of your geospatial data. For example, if you import all state boundaries from USAboundaries, you will be zoomed out wide enough to see AK and HI.
- If you’d rather have a tighter zoom, use fitBounds() to set a bounding box
Now your map will be perfect with no errors
To view your map in a bigger format than the viewer pane of Rstudio, click the Zoom or Show in New Window icons
You now have multiple options to save or publish your map
- Click Export and Save as Web Page
- Use saveWidget() from htmlwidgets package
  - These options will save .html file to you directory. You can share file and it opens in a web browser (Chrome best bet)
- Alternatively, publish map as web page! If you have set up an rpubs account you can do this in seconds by clicking Publish in viewer pane
  - Rpubs has 10-15 MB max size and is public
- Can also embed leaflet maps in R Markdown documents, shiny apps, etc.

The Map

This code is the same as is in the github repository. Just an annotated step by step explanation of what each code chunk is doing

# load packages -----------------------------------------------------------

## If a package is installed, it will be loaded. If any are not, the missing package(s) will be installed from CRAN and then loaded.

## First specify the packages of interest
packages <-  c('readxl', 'USAboundaries', 'sf', 'leaflet', 'leaflet.extras', 'htmlwidgets', 'tidyverse', 'tidylog')

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

# import data -------------------------------------------------------------

# we are pulling data from fwEco's github, so we need to direct R to that url so we don't need to type it everytime:
url <-  "https://raw.githubusercontent.com/FWeco/PAAFS_noStatic/master/"

# then we need to download each file necessary for mapping


# results dataframe - credit SRBC and USFWS NEFC - this was a spreadsheet, saved as rds for convenience
download.file(paste0(url, 'edna_map_df.rds'),'edna_map_df.rds', method="curl")
edna_map_df <- readRDS("edna_map_df.rds")
edna_map_df

# nhd Flowlines from shapefile - 7 components to each shapefile - download them all then import!
download.file(paste0(url, 'Flowlines.shp'), 'Flowlines.shp', method="curl")
download.file(paste0(url, 'Flowlines.cpg'), 'Flowlines.cpg', method="curl")
download.file(paste0(url, 'Flowlines.dbf'), 'Flowlines.dbf', method="curl")
download.file(paste0(url, 'Flowlines.prj'), 'Flowlines.prj', method="curl")
download.file(paste0(url, 'Flowlines.sbn'), 'Flowlines.sbn', method="curl")
download.file(paste0(url, 'Flowlines.sbx'), 'Flowlines.sbx', method="curl")
download.file(paste0(url, 'Flowlines.shx'), 'Flowlines.shx', method="curl")

nhd <- st_read('Flowlines.shp', stringsAsFactors = F) %>% # pipe!
  select(COMID, Resolution, GNIS_ID, GNIS_NAME, LENGTHKM, REACHCODE, FTYPE, StreamOrde, geometry) # select only cols we want

st_crs(nhd) # extract coordinate system (CS) information from an sf

nhd <- st_transform(nhd, "+proj=longlat +datum=WGS84") # transform projection to be WGS84 (same as other leaflet layers)

st_crs(nhd) # now in WGS84

nhd # view our 'attribute table'

## Simple feature collection with 303 features and 8 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -76.91974 ymin: 39.54369 xmax: -75.95318 ymax: 40.30795
## CRS:            +proj=longlat +datum=WGS84
## First 10 features:
##        COMID Resolution GNIS_ID         GNIS_NAME LENGTHKM      REACHCODE
## 1    4710026     Medium  591690 Susquehanna River    0.649 02050305000105
## 2    4710024     Medium  591690 Susquehanna River    0.157 02050305000105
## 3    4710020     Medium  591690 Susquehanna River    4.738 02050305000106
## 4    4710008     Medium  591690 Susquehanna River    1.946 02050305000107
## 5    4710012     Medium  591690 Susquehanna River    1.387 02050305000108
## 6    4710000     Medium  591690 Susquehanna River    4.250 02050305000109
## 7  932050170     Medium  591690 Susquehanna River    5.376 02050305000110
## 8    4710056     Medium  591690 Susquehanna River    0.277 02050305001444
## 9    4710032     Medium  591690 Susquehanna River    0.654 02050305007532
## 10   4710036     Medium  591690 Susquehanna River    4.327 02050305007534
##             FTYPE StreamOrde                       geometry
## 1  Artifical Path          7 LINESTRING (-76.79329 40.20...
## 2  Artifical Path          7 LINESTRING (-76.795 40.2021...
## 3  Artifical Path          7 LINESTRING (-76.84042 40.22...
## 4  Artifical Path          7 LINESTRING (-76.87124 40.23...
## 5  Artifical Path          7 LINESTRING (-76.85595 40.22...
## 6  Artifical Path          7 LINESTRING (-76.90595 40.26...
## 7  Artifical Path          7 LINESTRING (-76.91563 40.30...
## 8  Artifical Path          7 LINESTRING (-76.71948 40.12...
## 9  Artifical Path          7 LINESTRING (-76.78796 40.19...
## 10 Artifical Path          7 LINESTRING (-76.78162 40.19...

# dam shapefile 
download.file(paste0(url, 'eDNA_sampling_barriers.shp'), 'eDNA_sampling_barriers.shp', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.cpg'), 'eDNA_sampling_barriers.cpg', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.dbf'), 'eDNA_sampling_barriers.dbf', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.prj'), 'eDNA_sampling_barriers.prj', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.sbn'), 'eDNA_sampling_barriers.sbn', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.sbx'), 'eDNA_sampling_barriers.sbx', 
              method="curl")
download.file(paste0(url, 'eDNA_sampling_barriers.shx'), 'eDNA_sampling_barriers.shx', 
              method="curl")

barriers <- st_read('eDNA_sampling_barriers.shp', stringsAsFactors = F)

## Reading layer `eDNA_sampling_barriers' from data source `C:\R\trainings\2021-02-12_PAAFS_leaflet\fromGH\Rmd\eDNA_sampling_barriers.shp' using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 6 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -76.45289 ymin: 39.6146 xmax: -76.04288 ymax: 39.82653
## CRS:            4269

st_crs(barriers) # extract coordinate system (CS) information from an sf

## Coordinate Reference System:
##   User input: 4269 
##   wkt:
## GEOGCS["GCS_North_American_1983",
##     DATUM["North_American_Datum_1983",
##         SPHEROID["GRS_1980",6378137.0,298.257222101]],
##     PRIMEM["Greenwich",0.0],
##     UNIT["Degree",0.0174532925199433],
##     AUTHORITY["EPSG","4269"]]

barriers <- st_transform(barriers, "+proj=longlat +datum=WGS84") %>%  # transform projection to be WGS84 (same as other leaflet layers)
  mutate(Fish_Passa = ifelse(Fish_Passa == 'Y', 'Yes', 'No')) # make passage column more descriptive

## mutate: changed 7 values (100%) of 'Fish_Passa' (0 new NA)

st_crs(barriers) # now in WGS84

## Coordinate Reference System:
##   User input: +proj=longlat +datum=WGS84 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

First Map - Basic

# create legend based on values in presence column
sort(unique(edna_map_df$Presence))

## [1] "Absent"      "Not Sampled" "Present"

# create 
factpal_pres <- colorFactor(c('#91bfdb', '#bdbdbd', '#f03b20'), 
                            domain = edna_map_df$Presence, 
                            na.color = 'black')

# create basic map
leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap, group = "Basemap") %>% # basemap
  
  
  # add flowlines
  addPolylines(data = nhd,
               col = '#2b8cbe', weight = 2, opacity = 1,
               group = 'Impaired',
               label = ~ GNIS_NAME) %>% 
  
  
  # add barriers
  addCircleMarkers(
    data = barriers,
    lng = ~ Longitude, lat = ~ Latitude,
    col = ~ 'darkgrey', opacity = 0.8, 
    label = ~ Barrier) %>% 
  
  
  # add eDNA point results
  addCircleMarkers(
    data = edna_map_df,
    lng = ~ Longitude, lat = ~ Latitude,
    col = ~ factpal_pres(Presence), opacity = 1, 
    label = ~ Station,
    popup = ~ Presence) %>% 
  
  addLegend(pal = factpal_pres, values = edna_map_df$Presence, opacity = 0.8,
            title = 'Invasive Fish Presence')

Issues and Fixes

Flaw: little streams and big rivers all the same size
- Fix: vary nhd width by stream order variable
Flaw: can’t tell dams apart from eDNA samples
- Fix: images!
Flaw: can’t tell what species, or what year?
- Fix: Make separate layers for species and year
Flaw: Points that weren’t sampled are just as prominent as points that were
- Fix: Convert not sampled to NA, and set na color to transparent
Flaw: Dots/colors hard to see
- Fix: Fill circles with dark outline
Flaw: Would like to zoom in and see aerial imagery
- Fix: Add aerial basemap
Flaw: how to tell what county I’m in?
- Fix: add county boundaries
Flaw: how to orient myself once I’m way zoomed in?
- Fix: add a map reset button
Flaw: how to search for a specific town or set of coordinates?
- Fix: add a search box
Flaw: beginning extent too zoomed out
- Fix: add a bounding box to set extent of map zoom
Flaw: have no idea of distance between points/features
- Fix: add a scale bar

Advanced Map

# address fixes with data manipulation/include additional data files

# change 'Not Sampled' to NA
edna_map_df2 <- 
  edna_map_df %>% 
  mutate(Presence = ifelse(Presence == 'Not Sampled', NA_character_, Presence))
edna_map_df2

sort(unique(edna_map_df2$Presence))



# better legend with transparent NA's
# one for blue cats
factpal_bcat <- colorFactor(c('#91bfdb', '#fc8d59'), # color from HEX format - colorbrewer2.org
                             domain = edna_map_df2$Presence, 
                             na.color = 'transparent')

# one for snakeheads
factpal_nsh <- colorFactor(c('dodgerblue', 'brown4'), # or just colors from names
                            domain = edna_map_df2$Presence, 
                            na.color = 'transparent')


# get county boundaries from USAboundaries package
edna_co <- USAboundaries::us_counties(resolution = "high", states = c("PA", "MD")) %>%
  filter(name %in% c('York', 'Lancaster', 'Harford', 'Cecil'))

edna_co
st_crs(edna_co)
plot(edna_co$geometry)


# photo of a dam to use as barrier icon
damIcon <- makeIcon(
  iconUrl = paste0(url, 'dam.jpg'), # could use a web link
  iconWidth = 95, iconHeight = 38,
  iconAnchorX = 45, iconAnchorY = 37)

# purty map
edna_map <-
  
  leaflet() %>%
  
  addProviderTiles(providers$Esri.WorldTopoMap, group = "Topographic") %>% # basemap
  
  addProviderTiles(providers$Esri.WorldImagery, group = "Aerial") %>% # aerial map
  
  addResetMapButton() %>% # reset map extent on click
  
  addSearchOSM() %>% # search for towns/coords
  
  
  # add flowlines
  addPolylines(data = nhd, # sf data object with geometry field
               col = '#2b8cbe', # nice blue
               weight = ~ StreamOrde, # vary width of line by stream order!
               opacity = 1, # no transparency
               label = ~ GNIS_NAME, # show stream name on hover
               popup = ~ # make some fancy html popups that show on click
                 paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
                       '<br>',
                       '<b>', 'COMID: ', '</b>', COMID, "<br>",
                       '<b>', 'GNIS_ID:  ', '</b>', GNIS_ID, "<br>",
                       '<b>', 'GNIS_NAME:  ', '</b>', GNIS_NAME, "<br>",
                       '<b>', 'REACHCODE:  ', '</b>', REACHCODE, "<br>",
                       '<b>', 'Flowpath Type: ', '</b>', FTYPE, "<br>",
                       '<b>', 'Stream Order:  ', '</b>', StreamOrde, "<br>"),
               group = 'Hydrology') %>% # assign group so we can turn on/off
  
  
  # add Co boundary
  addPolygons(data = edna_co, # USAboundaries df
              color = "#737373", # light grey
              fill = F, # no fill inside polygon
              weight = 2) %>% # thickness of boundary
  
  
  # add barriers with jpg points!
  addMarkers(
    data = barriers, # not a sf object, just tabular
    lng = ~ Longitude, lat = ~ Latitude, # project lat/long
    icon = damIcon,
    popup = ~
      paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
            '<br>',
            '<b>', 'Barrier Name: ', '</b>', Barrier, "<br>",
            '<b>', 'Stream Name:  ', '</b>', Stream, "<br>",
            '<b>', 'Fish Passage?:  ', '</b>', Fish_Passa, "<br>",
            '<br>',
            '<b>', 'Latitude:  ', '</b>', Latitude, "<br>",
            '<b>', 'Longitude: ', '</b>', Longitude, "<br>",
            '<b>', 'Notes:  ', '</b>', Info, "<br>"),
    group = 'Dams') %>% 
  
  
  # points for NSH 2019
  addCircleMarkers(
    data = filter(edna_map_df2, Year == '2019' & Taxa == 'NSH'), # filter only 2019 NSH
    lng = ~ Longitude, lat = ~ Latitude, # project lat/long
    radius = 12, stroke = TRUE, # size of circle, and outer border
    col = '#252525', weight = 1, # color of border, thickness
    fill = TRUE, fillColor = ~ factpal_nsh(Presence), fillOpacity = 0.9, # fill with color from legend, 10% transparent
    label = ~ paste(Station, Year, Taxa, Presence, sep = '; '), # more advanced label on hover
    popup = ~
      paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
            '<br>',
            '<b>','<u>', paste(Year, Taxa, sep = ' - '), '</b>','</u>', "<br>",
            '<b>', Presence, '</b>', "<br>",
            '<br>',
            '<b>', 'Stream:  ', '</b>', Stream, "<br>",
            '<b>', 'Station:  ', '</b>', Station, "<br>",
            '<b>', 'Station-Filter:  ', '</b>', Station_Filter, "<br>",
            '<b>', 'Count: ', '</b>', round(Count, 1), "<br>",
            '<b>', 'Replicates:  ', '</b>', Replicates, "<br>", 
            '<b>', 'Volume Filtered:  ', '</b>', Volume, "<br>",
            '<br>',
            '<b>', 'Latitude: ', '</b>', Latitude, "<br>",
            '<b>', 'Longitude:  ', '</b>', Longitude, "<br>", 
            '<br>'),
    group = 'NSH 2019') %>% 
  
  
  # points for NSH 2020
  addCircleMarkers(
    data = filter(edna_map_df2, Year == '2020' & Taxa == 'NSH'),
    lng = ~ Longitude, lat = ~ Latitude,
    radius = 6, stroke = TRUE, col = '#252525', weight = 1, 
    fill = TRUE, fillColor = ~ factpal_nsh(Presence), fillOpacity = 0.9, 
    label = ~ paste(Station, Year, Taxa, Presence, sep = '; '),
    popup = ~
      paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
            '<br>',
            '<b>','<u>', paste(Year, Taxa, sep = ' - '), '</b>','</u>', "<br>",
            '<b>', Presence, '</b>', "<br>",
            '<br>',
            '<b>', 'Stream:  ', '</b>', Stream, "<br>",
            '<b>', 'Station:  ', '</b>', Station, "<br>",
            '<b>', 'Station-Filter:  ', '</b>', Station_Filter, "<br>",
            '<b>', 'Count: ', '</b>', round(Count, 1), "<br>",
            '<b>', 'Replicates:  ', '</b>', Replicates, "<br>", 
            '<b>', 'Volume Filtered:  ', '</b>', Volume, "<br>",
            '<br>',
            '<b>', 'Latitude: ', '</b>', Latitude, "<br>",
            '<b>', 'Longitude:  ', '</b>', Longitude, "<br>", 
            '<br>'),
    group = 'NSH 2020') %>% 
  
  
  # points for BCAT 2019
  addCircleMarkers(
    data = filter(edna_map_df2, Year == '2019' & Taxa == 'BCAT'),
    lng = ~ Longitude, lat = ~ Latitude,
    radius = 12, stroke = TRUE, col = '#252525', weight = 1, 
    fill = TRUE, fillColor = ~ factpal_bcat(Presence), fillOpacity = 0.6, 
    label = ~ paste(Station, Year, Taxa, Presence, sep = '; '),
    popup = ~
      paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
            '<br>',
            '<b>','<u>', paste(Year, Taxa, sep = ' - '), '</b>','</u>', "<br>",
            '<b>', Presence, '</b>', "<br>",
            '<br>',
            '<b>', 'Stream:  ', '</b>', Stream, "<br>",
            '<b>', 'Station:  ', '</b>', Station, "<br>",
            '<b>', 'Station-Filter:  ', '</b>', Station_Filter, "<br>",
            '<b>', 'Count: ', '</b>', round(Count, 1), "<br>",
            '<b>', 'Replicates:  ', '</b>', Replicates, "<br>", 
            '<b>', 'Volume Filtered:  ', '</b>', Volume, "<br>",
            '<br>',
            '<b>', 'Latitude: ', '</b>', Latitude, "<br>",
            '<b>', 'Longitude:  ', '</b>', Longitude, "<br>", 
            '<br>'),
    group = 'BCAT 2019') %>% 
  
  
  # points for BCAT 2020
  addCircleMarkers(
    data = filter(edna_map_df2, Year == '2020' & Taxa == 'BCAT'),
    lng = ~ Longitude, lat = ~ Latitude,
    radius = 6, stroke = TRUE, col = '#252525', weight = 1, 
    fill = TRUE, fillColor = ~ factpal_bcat(Presence), fillOpacity = 0.6, 
    label = ~ paste(Station, Year, Taxa, Presence, sep = '; '),
    popup = ~
      paste("<div class='leaflet-popup-scrolled' style='max-width:250px;max-height:250px'", 
            '<br>',
            '<b>','<u>', paste(Year, Taxa, sep = ' - '), '</b>','</u>', "<br>",
            '<b>', Presence, '</b>', "<br>",
            '<br>',
            '<b>', 'Stream:  ', '</b>', Stream, "<br>",
            '<b>', 'Station:  ', '</b>', Station, "<br>",
            '<b>', 'Station-Filter:  ', '</b>', Station_Filter, "<br>",
            '<b>', 'Count: ', '</b>', round(Count, 1), "<br>",
            '<b>', 'Replicates:  ', '</b>', Replicates, "<br>", 
            '<b>', 'Volume Filtered:  ', '</b>', Volume, "<br>",
            '<br>',
            '<b>', 'Latitude: ', '</b>', Latitude, "<br>",
            '<b>', 'Longitude:  ', '</b>', Longitude, "<br>", 
            '<br>'),
    group = 'BCAT 2020') %>% 
  
  
  # BCAT legend
  addLegend(position = "bottomright", pal = factpal_bcat, values = edna_map_df2$Presence, opacity = 0.9,
            title = 'Blue Catfish', group = c('BCAT 2019', 'BCAT 2020')) %>%

  # NSH legend
  addLegend(position = "bottomright", pal = factpal_nsh, values = edna_map_df2$Presence, opacity = 0.9,
            title = 'Northern Snakehead') %>% 
  
  
  # add scale bar
  addScaleBar(position = "bottomleft", options = scaleBarOptions(maxWidth = 200, metric = TRUE)) %>% 
  
  
  # Layers control
  addLayersControl(
    baseGroups = c("Topographic", "Aerial"),
    overlayGroups = c('NSH 2019', 'NSH 2020', 'BCAT 2019', 'BCAT 2020', 'Hydrology', 'Dams'), # groups that will overlay
    options = layersControlOptions(collapsed = FALSE)) %>% 
  
  hideGroup(c('NSH 2019', 'NSH 2020', 'BCAT 2019', 'BCAT 2020')) %>% # hide all groups at first, let user turn them on
  
  
  fitBounds(-76.5, 40.0, -75.9, 39.4) # bounding box so we are zoomed into good extent when map initializes

## filter: removed 204 rows (75%), 68 rows remaining
## filter: removed 204 rows (75%), 68 rows remaining
## filter: removed 204 rows (75%), 68 rows remaining
## filter: removed 204 rows (75%), 68 rows remaining

edna_map

She’s a beauty! Feel free to take this code and modify it for your use!

Hope this was helpful. If you care to read my R/science tweets follow @fwEco.

Pennsylvania AFS Workshop

Interactive Mapping in R: Don’t get caught up in the static

Matt Shank

February 12, 2021

The Biology

The Data Science

Using R as GIS

Tabular Data

Geospatial Data

Spatial Data Packages

The Pipe! `%>%`

Additional DS Skills

Introduction to leaflet

The Map

First Map - Basic

Issues and Fixes

Advanced Map

Pennsylvania AFS Workshop

Interactive Mapping in R: Don’t get caught up in the static

Matt Shank

February 12, 2021

The Biology

The Data Science

Using R as GIS

Tabular Data

Geospatial Data

Spatial Data Packages

The Pipe! %>%

Additional DS Skills

Introduction to leaflet

The Map

First Map - Basic

Issues and Fixes

Advanced Map

The Pipe! `%>%`