Invasive Species: Eucalyptus migration path

Loading the libraries

library(dismo)

## Loading required package: raster

## Loading required package: sp

library(maptools)

## Checking rgeos availability: TRUE

Using occ_search

We look for eucalyptus occurrence data in an online data library.

library(rgbif)
library(knitr)
eucalyptus_points <- occ_search(scientificName = "Eucalyptus", country = "MX", 
                                fields = c('name','country','countryCode','stateProvince',
                                           'year','decimalLatitude','decimalLongitude'),
                                limit = 300, return = 'data')

## Registered S3 method overwritten by 'crul':
##   method                 from
##   as.character.form_file httr

We make a simple map showing GPS points where occurrences of the species.

data(wrld_simpl)

xlim <- c(-129,-79)
ylim <- c(15,35)

plot(wrld_simpl, xlim = xlim, ylim = ylim)

points(eucalyptus_points$decimalLongitude, eucalyptus_points$decimalLatitude, col='red1')

We collect climate data from Mexico.

path <- file.path(system.file(package="dismo"), 'ex')
files <- list.files(path, pattern='grd$', full.names=TRUE )
files

## [1] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio1.grd" 
## [2] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio12.grd"
## [3] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio16.grd"
## [4] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio17.grd"
## [5] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio5.grd" 
## [6] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio6.grd" 
## [7] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio7.grd" 
## [8] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/bio8.grd" 
## [9] "C:/Users/Humberto Medina/OneDrive/Documentos/R/win-library/3.6/dismo/ex/biome.grd"

predictors <- stack(files)
predictors

## class      : RasterStack 
## dimensions : 192, 186, 35712, 9  (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5  (x, y)
## extent     : -125, -32, -56, 40  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names      : bio1, bio12, bio16, bio17, bio5, bio6, bio7, bio8, biome 
## min values :  -23,     0,     0,     0,   61, -212,   60,  -66,     1 
## max values :  289,  7682,  2458,  1496,  422,  242,  461,  323,    14

extent_of_mexico_map <- extent(-129, -79, -15, 35)
predictors_cropped_to_mexico <- crop(predictors, extent_of_mexico_map)

set.seed(0)
group <- kfold(eucalyptus_points, 5)
pres_train_your_species <- eucalyptus_points[group != 1, ]
pres_train_your_species <- as.data.frame(pres_train_your_species[,1:2])

pres_test_your_species <- eucalyptus_points[group == 1, ]
pres_test_your_species <- as.data.frame(pres_test_your_species[,1:2])

pred_nf <- dropLayer(predictors_cropped_to_mexico, 'biome')
backg <- randomPoints(pred_nf, n=1000, ext=extent_of_mexico_map, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

We make the MAxEnt dysfunction model

jar <- paste(system.file(package="dismo"), "/java/maxent.jar", sep='')
xm <- maxent(predictors_cropped_to_mexico, pres_train_your_species, factors='biome')

## Warning in .local(x, p, ...): 1 (1.64%) of the presence points have NA
## predictor values

## Loading required namespace: rJava

plot(xm)

We build the map.

e <- evaluate(pres_test_your_species, backg_test, xm, predictors_cropped_to_mexico)
e

## class          : ModelEvaluation 
## n presences    : 52 
## n absences     : 199 
## AUC            : 0.9527445 
## cor            : 0.7425377 
## max TPR+TNR at : 0.2549756

px <- predict(predictors_cropped_to_mexico, xm, ext=extent_of_mexico_map, progress='')
par(mfrow=c(1,2))
plot(px, main='## Eucalyptus ##')

Appearances by state:

states_factor <- as.factor(eucalyptus_points$stateProvince)
counted_states <- plyr::count(states_factor)
plot(counted_states$x, counted_states$freq)