Here we provide the code for the article: Elevation as occupancy determinant of the little red brocket deer (Mazama rufina) in the Central Andes of Colombia.
library(raster)
library(sf)
library(unmarked)
library(tidyverse)
library(lubridate)
library(readxl)
library(mapview)
library (tmap) # nice plane maps
library (tmaptools) # more maps
library (osmdata) # read osm
library (OpenStreetMap) # osm maps
library (grid) # mix maps
library (GADMTools) # subset GADM
library (rasterVis)# Quindio_data <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Vanessa/data/matriz_v3.xlsx")
Quindio_data <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/matriz_v4.xlsx")
Quindio_data$Longitud <- Quindio_data$Longitud * (-1)
Risaralda_data <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/matriz_v2.xlsx", sheet = "Ucumari"))
# covaribles
# Covs <- read.delim("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/Covariables.csv")
#
# Covs_q <- data.frame(scale(Covs[,5:12])) # estandarizadas
# Risaralda_data <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/matriz_v2.xlsx", sheet = "Ucumari"))
Risaralda_data [,4] <- as.vector(as.numeric(Risaralda_data [,4]))
Risaralda_data [,5] <- as.vector(as.numeric(Risaralda_data [,5]))
Risaralda_data [,17] <- as.vector(as.numeric(Risaralda_data [,17]))
Risaralda_data$binomial <- paste(Risaralda_data$Genus, Risaralda_data$Species , sep = "_")
Risaralda_data$camera_trap_start_date <- ymd(Risaralda_data$camera_trap_start_date)
Risaralda_data$camera_trap_end_date <- ymd(Risaralda_data$camera_trap_end_date)Risaralda_sf <- st_as_sf(Risaralda_data, coords = c("Longitude", "Latitude"), crs = "+proj=longlat +datum=WGS84 +no_defs")
cams_plot <- mapview(Risaralda_sf, zcol = "Photo Type",
color = c("black","red"),
legend = FALSE,
map.types = "Esri.WorldImagery")
Quindio_sf <- st_as_sf(Quindio_data, coords = c("Longitud", "Latitud"), crs = "+proj=longlat +datum=WGS84 +no_defs")
Risaralda_data$camara <- Risaralda_data$Camera_Trap
Risaralda_data$Latitud <- Risaralda_data$Latitude
Risaralda_data$Longitud <- Risaralda_data$Longitude
# cams_plot
### Risaralda
cam_location_R <- unique(Risaralda_data[,c('Latitud','Longitud','camara')])
#as.data.frame(unique(cbind(Risaralda_data$Longitude,
#Risaralda_data$Latitude)))
#### Quindio #### las Camaras se sobrelapan
cam_location_Q <- unique(Quindio_data[,c('Latitud','Longitud','camara')])
cams_loc_QR <- rbind(cam_location_R, cam_location_Q)
cams_loc_QR_sf <- st_as_sf(cams_loc_QR, coords = c("Longitud", "Latitud"), crs = "+proj=longlat +datum=WGS84 +no_defs")
cams_sp_R <- st_as_sf(cam_location_R, coords = c("Longitud", "Latitud"), crs = "+proj=longlat +datum=WGS84 +no_defs")
centroid <- c(mean(Risaralda_data$Longitud), mean(Risaralda_data$Latitud))
cams_sp_Q <- as(Quindio_sf, "Spatial")
##################################
# set extent get SRTM
##################################
clip_window <- extent(-75.60 , -75.39, 4.59, 4.81)
bb <- c(-75.60, 4.59, -75.39, 4.81)
srtm <- raster::getData('SRTM', lon=centroid[1], lat=centroid[2])
# plot(clip_window, border = "blue")
# plot(cams_sp, add=TRUE)
# crop the raster using the vector extent
srtm_crop <- crop(srtm, clip_window)
# elevation.crop
elevation <- srtm_crop
slope<-terrain(srtm_crop, opt="slope", unit='degrees', neighbors=8)
aspect<-terrain(srtm_crop, opt="aspect", unit='degrees', neighbors=8)
roughness <- terrain(srtm_crop, opt = c("roughness"))
cov.stack<-stack(elevation, slope, aspect, roughness)
names(cov.stack) <- c("elevation", "slope", "aspect", "roughness" )
# extract altitudes
cam_covs <- raster::extract(cov.stack, cams_loc_QR_sf)
full_covs <- as.data.frame(cam_covs) # convert to Data frame
full_covs_1 <- scale(full_covs)
full_covs_s <- as.data.frame(full_covs_1)
full_covs_s$camara <- cams_loc_QR$camara
# plot(srtm_crop)
# plot(cams_sp, add=TRUE)
# view
pal = mapviewPalette("mapviewTopoColors")
# Risaralda_sf %>% group_by(Camera_Trap) %>%
# summarise(mean = mean("Number of Animals", na.omit=T), fotos = n()) %>%
# filter(fotos > 1) %>%
# mapview (cex = "fotos", zcol = "fotos", legend = TRUE) +
# mapview (srtm_crop, col.regions = pal(400), at = seq(1000, 5000, 10), alpha=0.7) +
# mapview(Risaralda_sf, cex = 1) + mapview(st_jitter(Quindio_sf, factor = 0.03), alpha = 0, cex = 4, zcol = "camara")
#
# get fondo de osm
andes_osm1 <- read_osm(bb, zoom = NULL, type="stamen-terrain", mergeTiles = TRUE) # type puede ser tambien bing, osm # type puede ser tambien bing, osm
# get fondo de osm
# col_osm1 <- read_osm(bb2, type="stamen-toner") # type puede ser tambien bing, osm
# col <- raster::getData('GADM', country='COL', level=1)
# col.sf <- st_as_sf(col, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
colombia <- gadm_sf_loadCountries("COL", level=1, basefile="./")
collimit <- gadm_sf_loadCountries("COL", level=0, basefile="./")
deptos <- gadm_subset(colombia, regions=c("Risaralda", "QuindÃo"))
depto_window <- qtm(andes_osm1) +
tm_shape(cams_loc_QR_sf) +
tm_dots(col = "red", size = 0.2,
shape = 16, title = "Sampling point", legend.show = TRUE,
legend.is.portrait = TRUE,
legend.z = NA) +
tm_layout(scale = .9) +
# legend.position = c(.78,.72),
# legend.outside.size = 0.1,
# legend.title.size = 1.6,
# legend.height = 0.9,
# legend.width = 1.5,
# legend.text.size = 1.2) +
# legend.hist.size = 0.5) +
tm_legend(position = c("left", "bottom"), frame = TRUE,
bg.color="white") +
tm_layout(frame=F) + tm_scale_bar() + tm_compass(position = c(.75, .82), color.light = "grey90")
dep_map <- tm_shape(deptos$sf) + tm_polygons() +
tm_shape(Risaralda_sf) + tm_symbols(shape = 0, col = "red", size = 0.25)
col_map <- tm_shape(collimit$sf) + tm_polygons() + tm_shape(deptos$sf) + tm_polygons()
##### print ambos
depto_window
print(dep_map, vp = viewport(0.73, 0.40, width = 0.25, height = 0.25))
print(col_map, vp = viewport(0.73, 0.65, width = 0.25, height = 0.25))# ggplot() + geom_sf(data = Risaralda_sf) + geom_sf(data = Risaralda_sf,
# aes(color = Genus ),
# alpha = 0.7)
# show.legend = "point")
# Risaralda_sf %>% group_by(Camera_Trap) %>%
# summarise(mean = mean("Number of Animals", na.omit=T), n = n())%>%
# mapview(cex = n)
plot(cov.stack)# source("C:/Users/diego.lizcano/Box Sync/CodigoR/Vanessa/R/team.R")
source("D:/BoxFiles/Box Sync/CodigoR/Vanessa/R/team.R")
mat.per.sp<-f.matrix.creator2(data = Risaralda_data, year = c(2017, 2016))
sp.names<-names(mat.per.sp) # species names
# counting how many (total) records per species by all days
cont.per.sp<-data.frame(row.names = sp.names)
for (i in 1:length(mat.per.sp)){
cont.per.sp[i,1]<-sum(apply(as.data.frame(mat.per.sp [[i]]),FUN=sum,na.rm=T, MARGIN = 1))
}
cont.per.sp$especie<-rownames(cont.per.sp)
colnames(cont.per.sp)<-c("Numero_de_fotos","especie")
# nombre de la matriz de cada especie
# names(mat.per.sp)
print(as.data.frame(arrange(df = cont.per.sp, desc(Numero_de_fotos))))## Numero_de_fotos especie
## 1 53 NA_NA
## 2 16 Mazama_rufina
Quindio_data$Elev <- cam_location_Q$Elev_Q # pega elevation
y_Q <- as.data.frame(Quindio_data[,4:111])
y_R <- as.data.frame(mat.per.sp[[1]])
# Covariables
#Correr el modelo de ocupacion
umf_Q<- unmarkedFrameOccu(y= y_Q, siteCovs = NULL) # Quindio
umf_R<- unmarkedFrameOccu(y= y_R, siteCovs = NULL) # Risaralda
#######Graficar umf
plot(umf_Q, main="Quindio")## unmarkedFrame Object
##
## 28 sites
## Maximum number of observations per site: 108
## Mean number of observations per site: 34.25
## Sites with at least one detection: 11
##
## Tabulation of y observations:
## 0 1 <NA>
## 938 21 2065
## unmarkedFrame Object
##
## 59 sites
## Maximum number of observations per site: 96
## Mean number of observations per site: 45.32
## Sites with at least one detection: 7
##
## Tabulation of y observations:
## 0 1 <NA>
## 2658 16 2990
#Construir modelos
# nombre camara en quindio
row.names(y_Q) <- Quindio_data$camara
# nuevo Risaralda con mas NAs en columnas
y_Rn <- cbind(y_R, as.data.frame(matrix(NA,59,12))) # Nuevo Risaralda con 12 NAs
#names(y_Rn) <- NULL
#names(y_Q) <- NULL
# Loop to convert logical to numeric
for (i in 1:length(y_Rn)){
y_Rn[,i] <- as.numeric(y_Rn[,i])
}
names(y_Rn) <- c(1:108)
# Loop to convert logical to numeric
for (i in 1:length(y_Q)){
y_Q[,i] <- as.numeric(y_Q[,i])
}
names(y_Q) <- c(1:108)
# unir Risaralda y Quindio
y_full_d <- rbind (y_Rn, y_Q)
y_full_d$camara <- row.names(y_full_d)
# sort both per camera
y_full <- y_full_d[order(y_full_d$camara),]
full_covs_s2 <- full_covs_s[order(full_covs_s$camara),]
# Make unmarked frame
umf_y_full<- unmarkedFrameOccu(y= y_full[,1:108])
siteCovs(umf_y_full) <- full_covs_s # data.frame(Elev=full_covs$Elev) # Full
#######Graficar umf
plot(umf_y_full)# build models
mf0<-occu(~1 ~ 1, umf_y_full)
mf1<-occu(~1 ~ elevation, umf_y_full)
mf2<-occu(~1 ~ elevation +I(elevation^2), umf_y_full)
mf3<-occu(~1 ~ slope, umf_y_full)
mf4<-occu(~1 ~ aspect, umf_y_full)
mf5<-occu(~1 ~ roughness, umf_y_full, starts = c(0.6, -3, 0))
mf6<-occu(~elevation +I(elevation^2) ~ elevation +I(elevation^2), umf_y_full)
mf7<-occu(~roughness ~ elevation +I(elevation^2), umf_y_full)
mf8<-occu(~slope ~ elevation +I(elevation^2), umf_y_full)
# fit list
fms1<-fitList("p(.) Ocu(.)"=mf0,
"p(.) Ocu(elev)"=mf1,
"p(.) Ocu(elev^2)"=mf2,
"p(.) Ocu(slope)"=mf3,
"p(.) Ocu(aspect)"=mf4,
"p(.) Ocu(roughness)"=mf5,
"p(elev^2) Ocu(elev^2)"=mf6,
"p(roughness) Ocu(elev^2)"=mf7,
"p(slope) Ocu(elev^2)"=mf8
)
modSel(fms1)## nPars AIC delta AICwt cumltvWt
## p(roughness) Ocu(elev^2) 5 359.00 0.00 8.4e-01 0.84
## p(slope) Ocu(elev^2) 5 362.46 3.46 1.5e-01 0.99
## p(.) Ocu(elev^2) 4 368.33 9.33 7.9e-03 0.99
## p(.) Ocu(elev) 3 369.73 10.73 3.9e-03 1.00
## p(elev^2) Ocu(elev^2) 6 372.25 13.25 1.1e-03 1.00
## p(.) Ocu(aspect) 3 374.04 15.03 4.6e-04 1.00
## p(.) Ocu(roughness) 3 378.38 19.37 5.2e-05 1.00
## p(.) Ocu(.) 2 417.05 58.05 2.1e-13 1.00
## p(.) Ocu(slope) 3 419.05 60.05 7.7e-14 1.00
## t0 = 36.31634
newdat_range<-data.frame(elevation=seq(min(full_covs_s$elevation),
max(full_covs_s$elevation),length=100),
roughness=seq(min(full_covs_s$roughness),
max(full_covs_s$roughness), length=100))
## plot Detection en escala estandarizada
pred_det <-predict(mf7, type="det", newdata=newdat_range, appendData=TRUE)
plot(Predicted~roughness, pred_det,type="l",col="blue",
xlab="roughness",
ylab="Probabilidad de detección")
lines(lower~roughness, pred_det,type="l",col=gray(0.5))
lines(upper~roughness, pred_det,type="l",col=gray(0.5))### plot occupancy en escala estandarizada
pred_psi <-predict(mf7, type="state", newdata=newdat_range, appendData=TRUE)
plot(Predicted~elevation, pred_psi,type="l",col="blue", ylim=c(0,0.95),
xlab="Altitud",
ylab="Probabilidad de ocupación")
lines(lower~elevation, pred_psi,type="l",col=gray(0.5))
lines(upper~elevation, pred_psi,type="l",col=gray(0.5))## plot Detection en escala original
pred_det <-predict(mf7, type="det", newdata=newdat_range, appendData=TRUE)
plot(Predicted~roughness, pred_det,type="l",col="blue",
xlab="Roughness",
ylab="Detection Probability",
xaxt="n")
xticks <- c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3) # -1:2
xlabs <- xticks*sd(full_covs$roughness) + mean(full_covs$roughness) #Use the mean and sd of original value to change label name
axis(1, at=xticks, labels=round(xlabs, 1))
lines(lower~roughness, pred_det,type="l",col=gray(0.5))
lines(upper~roughness, pred_det,type="l",col=gray(0.5))### Plot occupancy en escala original
plot(Predicted ~ elevation, pred_psi, type="l", ylim=c(0,1), col="blue",
xlab="Elevation",
ylab="Occupancy Probability",
xaxt="n")
xticks <- c(-1, -0.5, 0, 0.5, 1, 1.5, 2) # -1:2
xlabs <- xticks*sd(full_covs$elevation) + mean(full_covs$elevation) #Use the mean and sd of original value to change label name
axis(1, at=xticks, labels=round(xlabs, 1))
lines(lower ~ elevation, pred_psi, type="l", col=gray(0.5))
lines(upper ~ elevation, pred_psi, type="l", col=gray(0.5))La altitud es buena predictora de la ocupación. la ocupación aumenta con la altitud hasta los 3100 metros, luego de esta altitud, la ocupación disminuye.
library(RColorBrewer)
srtm_crop_s <- stack(scale(elevation),
scale(roughness)) # scale altitud
names(srtm_crop_s) <- c("elevation", "roughness")
crs(srtm_crop_s) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
pred_p_s <-predict(mf7, type="det", newdata=srtm_crop_s) ## doing row 1000 of 66528
## doing row 2000 of 66528
## doing row 3000 of 66528
## doing row 4000 of 66528
## doing row 5000 of 66528
## doing row 6000 of 66528
## doing row 7000 of 66528
## doing row 8000 of 66528
## doing row 9000 of 66528
## doing row 10000 of 66528
## doing row 11000 of 66528
## doing row 12000 of 66528
## doing row 13000 of 66528
## doing row 14000 of 66528
## doing row 15000 of 66528
## doing row 16000 of 66528
## doing row 17000 of 66528
## doing row 18000 of 66528
## doing row 19000 of 66528
## doing row 20000 of 66528
## doing row 21000 of 66528
## doing row 22000 of 66528
## doing row 23000 of 66528
## doing row 24000 of 66528
## doing row 25000 of 66528
## doing row 26000 of 66528
## doing row 27000 of 66528
## doing row 28000 of 66528
## doing row 29000 of 66528
## doing row 30000 of 66528
## doing row 31000 of 66528
## doing row 32000 of 66528
## doing row 33000 of 66528
## doing row 34000 of 66528
## doing row 35000 of 66528
## doing row 36000 of 66528
## doing row 37000 of 66528
## doing row 38000 of 66528
## doing row 39000 of 66528
## doing row 40000 of 66528
## doing row 41000 of 66528
## doing row 42000 of 66528
## doing row 43000 of 66528
## doing row 44000 of 66528
## doing row 45000 of 66528
## doing row 46000 of 66528
## doing row 47000 of 66528
## doing row 48000 of 66528
## doing row 49000 of 66528
## doing row 50000 of 66528
## doing row 51000 of 66528
## doing row 52000 of 66528
## doing row 53000 of 66528
## doing row 54000 of 66528
## doing row 55000 of 66528
## doing row 56000 of 66528
## doing row 57000 of 66528
## doing row 58000 of 66528
## doing row 59000 of 66528
## doing row 60000 of 66528
## doing row 61000 of 66528
## doing row 62000 of 66528
## doing row 63000 of 66528
## doing row 64000 of 66528
## doing row 65000 of 66528
## doing row 66000 of 66528
pred_p_r <- pred_p_s * sd(full_covs$roughness) + mean(full_covs$roughness)
crs(pred_p_r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
clr <- colorRampPalette(brewer.pal(9, "YlGn"))
mapview (pred_p_r[[1]], col.regions = clr, legend = TRUE, alpha=0.7)# plot(pred_psi_s[[1]], main="Detection", col = topo.colors(20))
levelplot(pred_p_s[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Detection")library(RColorBrewer)
# srtm_crop_s <- stack(scale(elevation),
# scale(roughness)) # scale altitud
# names(srtm_crop_s) <- c("elevation", "roughness")
# crs(srtm_crop_s) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
pred_psi_s <-predict(mf7, type="state", newdata=srtm_crop_s) ## doing row 1000 of 66528
## doing row 2000 of 66528
## doing row 3000 of 66528
## doing row 4000 of 66528
## doing row 5000 of 66528
## doing row 6000 of 66528
## doing row 7000 of 66528
## doing row 8000 of 66528
## doing row 9000 of 66528
## doing row 10000 of 66528
## doing row 11000 of 66528
## doing row 12000 of 66528
## doing row 13000 of 66528
## doing row 14000 of 66528
## doing row 15000 of 66528
## doing row 16000 of 66528
## doing row 17000 of 66528
## doing row 18000 of 66528
## doing row 19000 of 66528
## doing row 20000 of 66528
## doing row 21000 of 66528
## doing row 22000 of 66528
## doing row 23000 of 66528
## doing row 24000 of 66528
## doing row 25000 of 66528
## doing row 26000 of 66528
## doing row 27000 of 66528
## doing row 28000 of 66528
## doing row 29000 of 66528
## doing row 30000 of 66528
## doing row 31000 of 66528
## doing row 32000 of 66528
## doing row 33000 of 66528
## doing row 34000 of 66528
## doing row 35000 of 66528
## doing row 36000 of 66528
## doing row 37000 of 66528
## doing row 38000 of 66528
## doing row 39000 of 66528
## doing row 40000 of 66528
## doing row 41000 of 66528
## doing row 42000 of 66528
## doing row 43000 of 66528
## doing row 44000 of 66528
## doing row 45000 of 66528
## doing row 46000 of 66528
## doing row 47000 of 66528
## doing row 48000 of 66528
## doing row 49000 of 66528
## doing row 50000 of 66528
## doing row 51000 of 66528
## doing row 52000 of 66528
## doing row 53000 of 66528
## doing row 54000 of 66528
## doing row 55000 of 66528
## doing row 56000 of 66528
## doing row 57000 of 66528
## doing row 58000 of 66528
## doing row 59000 of 66528
## doing row 60000 of 66528
## doing row 61000 of 66528
## doing row 62000 of 66528
## doing row 63000 of 66528
## doing row 64000 of 66528
## doing row 65000 of 66528
## doing row 66000 of 66528
pred_psi_r <- pred_psi_s * sd(full_covs$elevation) + mean(full_covs$elevation)
crs(pred_psi_r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
clr <- colorRampPalette(brewer.pal(9, "YlGn"))
mapview (pred_psi_r[[1]], col.regions = clr, legend = TRUE, alpha=0.7)# plot(pred_psi_s[[1]], main="Occupancy")
levelplot(pred_psi_r[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] grid parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 plyr_1.8.5 rasterVis_0.47
## [4] latticeExtra_0.6-29 GADMTools_3.7-2 rgdal_1.4-8
## [7] classInt_0.4-2 OpenStreetMap_0.3.4 osmdata_0.1.3
## [10] tmaptools_2.0-2 tmap_2.3-2 mapview_2.7.0
## [13] readxl_1.3.1 lubridate_1.7.4 forcats_0.4.0
## [16] stringr_1.4.0 dplyr_0.8.4 purrr_0.3.3
## [19] readr_1.3.1 tidyr_1.0.2 tibble_2.1.3
## [22] ggplot2_3.2.1 tidyverse_1.3.0 unmarked_0.13-1
## [25] reshape2_1.4.3 Rcpp_1.0.3 lattice_0.20-38
## [28] sf_0.8-1 raster_3.0-12 sp_1.4-0
##
## loaded via a namespace (and not attached):
## [1] leafem_0.0.1 colorspace_1.4-1 rjson_0.2.20
## [4] class_7.3-15 leaflet_2.0.3 satellite_1.0.2
## [7] base64enc_0.1-3 fs_1.3.1 dichromat_2.0-0
## [10] rstudioapi_0.11 farver_2.0.3 hexbin_1.28.1
## [13] fansi_0.4.1 xml2_1.2.2 codetools_0.2-16
## [16] rosm_0.2.5 knitr_1.28 jsonlite_1.6.1
## [19] rJava_0.9-11 broom_0.5.4 dbplyr_1.4.2
## [22] png_0.1-7 rgeos_0.5-2 shiny_1.4.0
## [25] compiler_3.6.1 httr_1.4.1 backports_1.1.5
## [28] assertthat_0.2.1 Matrix_1.2-18 fastmap_1.0.1
## [31] lazyeval_0.2.2 cli_2.0.1 later_1.0.0
## [34] leaflet.providers_1.9.0 htmltools_0.4.0 tools_3.6.1
## [37] ggmap_3.0.0 gtable_0.3.0 glue_1.3.1
## [40] cellranger_1.1.0 vctrs_0.2.2 svglite_1.2.3
## [43] nlme_3.1-144 leafsync_0.1.0 crosstalk_1.0.0
## [46] lwgeom_0.2-1 xfun_0.12 rvest_0.3.5
## [49] mime_0.9 lifecycle_0.1.0 XML_3.99-0.3
## [52] zoo_1.8-7 MASS_7.3-51.5 scales_1.1.0
## [55] hms_0.5.3 promises_1.1.0 yaml_2.2.1
## [58] curl_4.3 gridExtra_2.3 gdtools_0.2.1
## [61] ggspatial_1.0.3 stringi_1.4.6 maptools_0.9-9
## [64] leafpop_0.0.5 e1071_1.7-3 systemfonts_0.1.1
## [67] bitops_1.0-6 RgoogleMaps_1.4.5.3 rlang_0.4.4
## [70] pkgconfig_2.0.3 evaluate_0.14 htmlwidgets_1.5.1
## [73] tidyselect_1.0.0 magrittr_1.5 R6_2.4.1
## [76] generics_0.0.2 DBI_1.1.0 pillar_1.4.3
## [79] haven_2.2.0 foreign_0.8-75 withr_2.1.2
## [82] units_0.6-5 modelr_0.1.5 crayon_1.3.4
## [85] uuid_0.1-2 KernSmooth_2.23-16 rmarkdown_2.1
## [88] jpeg_0.1-8.1 reprex_0.3.0 digest_0.6.24
## [91] webshot_0.5.2 xtable_1.8-4 httpuv_1.5.2
## [94] brew_1.0-6 stats4_3.6.1 munsell_0.5.0
## [97] viridisLite_0.3.0