Interactive map with leaflet and R: plotting a large quantity of sites
We will see here how to plot interactively a huuuuuge amount of site on a map with very few lines of code. I am not an expert in GIS, but people have developed nice packages that are doing the job for you. At the end of this tutorial, you will be able to do these maps (and much better ones) by just copy/pasting my code and change the input data.
The map was published here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8507310/bin/13071_2021_4950_MOESM1_ESM.html
In the paper from Roland Bamou et al. An update on the mosquito fauna and mosquito-borne diseases distribution in Cameroon, 2021, Parasites and vectors.
Let’s go.
You will need to install and load these packages:
-leaflet
-plyr
-reshape2
-RColorBrewer
-stringr
Load data
dat <- read.table("Database_ano.txt", h=T )
head(dat)## Country GAUL_Admin1 GAUL_Admin2 Lat Long YeStart YeEnd
## 1 Cameroon Adamaoua Djerem 6.500000 12.70000 1998 2007
## 2 Cameroon Adamaoua Djerem 6.509783 12.47790 1998 2007
## 3 Cameroon Adamaoua Djerem 6.569880 12.69280 1998 2007
## 4 Cameroon Adamaoua Djerem 6.751020 13.11061 1998 2007
## 5 Cameroon Adamaoua Djerem 6.366667 12.75000 1998 2007
## 6 Cameroon Adamaoua Djerem 6.550000 12.68333 1998 2007
## An_gambiae_complex An_gambiae_ss SS_M_Form SS_S_Form An_arabiensis An._melas
## 1 1 1 0 0 0 0
## 2 1 1 0 1 1 0
## 3 1 1 0 1 1 0
## 4 1 1 1 1 1 0
## 5 1 1 0 0 0 0
## 6 1 1 0 0 0 0
## An._merus An_bwambae An_funestus_s.l An_funestus_s.s. An_rivulorum An_leesoni
## 1 0 0 1 1 0 0
## 2 0 0 1 1 0 0
## 3 0 0 1 1 0 0
## 4 0 0 1 1 0 0
## 5 0 0 1 1 0 0
## 6 0 0 1 1 0 0
## An_parensis An_vaneedeni An_nili_s.l An_moucheti_s.l An_pharoensis
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
## An_hancocki An_mascarensis An_marshalli An_squamous An_wellcomei An_rufipes
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## An_coustani_s.l An_ziemanni An_paludis Adults_Larvae Sampling_Methods
## 1 0 0 0 Adults PSC
## 2 0 0 0 Adults PSC
## 3 0 0 0 Adults PSC
## 4 0 0 0 Adults PSC
## 5 0 0 0 Adults PSC
## 6 0 0 0 Adults PSC
## Species_Identification
## 1 M_PCR
## 2 M_PCR
## 3 M_PCR
## 4 M_PCR
## 5 M_PCR
## 6 M_PCRPrepare the data
I melt the data from wide to long data.frame format, discard some species that I don’t want to plot and make some text edit.
data_anopheles <- melt(dat[,-c(8,34,35,36)], c(1:7))
data_anopheles <- na.omit(data_anopheles)
data_anopheles$Lat <- as.numeric(as.character(data_anopheles$Lat))
data_anopheles$Long <- as.numeric(as.character(data_anopheles$Long))
ddply(data_anopheles, .(variable), summarise, N=sum(value) )## variable N
## 1 An_gambiae_ss 456
## 2 SS_M_Form 137
## 3 SS_S_Form 318
## 4 An_arabiensis 227
## 5 An._melas 11
## 6 An._merus 0
## 7 An_bwambae 0
## 8 An_funestus_s.l 439
## 9 An_funestus_s.s. 232
## 10 An_rivulorum 11
## 11 An_leesoni 13
## 12 An_parensis 0
## 13 An_vaneedeni 0
## 14 An_nili_s.l 141
## 15 An_moucheti_s.l 153
## 16 An_pharoensis 54
## 17 An_hancocki 59
## 18 An_mascarensis 0
## 19 An_marshalli 21
## 20 An_squamous 26
## 21 An_wellcomei 23
## 22 An_rufipes 41
## 23 An_coustani_s.l 132
## 24 An_ziemanni 44
## 25 An_paludis 66data_anophelesb <- subset(data_anopheles, variable!= "An._merus" & variable!= "An_bwambae" & variable!= "An_parensis"
& variable!= "An_vaneedeni"
& variable!= "An_mascarensis"
& variable!= "An_gambiae_ss")
#ddply(data_anophelesb, .(variable), summarise, N=sum(value) )
#Replace text for a better visual
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "SS_M_Form", "An. gambiae ss M form")
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "SS_S_Form", "An. gambiae ss S form")
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "_", ". ")
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "An.. melas", "An. melas")
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "_s.l", " sl")
data_anophelesb$variable <- str_replace(data_anophelesb$variable, "_s.s.", " ss")Create your Bounding Box and visualize with a Leaflet interactive map.
GPS coordinates are in Decimal Degrees. If your coordinates are in DMS format (Degrees Minutes Seconds), you can easily convert them via dedicated websites, such as: https://www.latlong.net/lat-long-dms.html (this one display the location on a map).
# Define bounding box with longitude/latitude coordinates
bbox <- list(
p1 = list(long = 8, lat = 1),
p2 = list(long = 17, lat = 14)
)Is the bounding box ok? Change the GPS coordinate if needed.
Define a clolor palette
I assign a color to each species.
n <- 20
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
pie(rep(1,n), col=sample(col_vector, n))
data_anophelesb$variable <- factor(data_anophelesb$variable)
Col <- colorFactor(palette = col_vector, data_anophelesb$variable)I now only select sites where a species was observed
data_anophelesb <- subset(data_anophelesb, value==1 )Plot the map with Leaflet
Define your group factor
Here the grouping factor is the Anopheles species.
groups = as.character(unique(data_anophelesb$variable))Make the map. The name of the place will appear when you click on a point on the map (popup), sample date will appear when you will place the mouse over the point (label).
map = leaflet() %>%
addTiles() %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
addRectangles(
lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
fillColor = "transparent" ) %>%
fitBounds(
lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
lng2 = bbox$p2$long, lat2 = bbox$p2$lat )
for(g in groups){
d = data_anophelesb[data_anophelesb$variable == g, ]
map = map %>% addCircleMarkers(data = d, ~Long, ~Lat, popup = ~as.character(GAUL_Admin2),color = ~Col(variable), label = ~as.character(paste(d$YeStart,d$YeEnd, sep="-")), radius=5, stroke = FALSE, fillOpacity = 0.9, group = g)
}Add layer controls and a minimap and display the results:
map %>% addLayersControl(baseGroups = c("OSM (default)", "Toner Lite"), overlayGroups = groups, options = layersControlOptions(collapsed = FALSE)) %>% addMiniMap(tiles = providers$Esri.OceanBasemap, width = 120, height=80)You can now save a stand alone htlm file from Rstudio with export/Save a web page.
NOTE: The layer controls are not rendering here. But it should appear OK on your computer. Here is the final map: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8507310/bin/13071_2021_4950_MOESM1_ESM.html
Albin Fontaine
Researcher
My main research interests include vector-borne diseases surveillance, interaction dynamics between mosquito-borne pathogens and their hosts, and quantitative genetics applied to mosquito-viruses interactions.