Pachystachya - distribution modeling
Nguyễn Lê Xuân Bách
April 6, 2017
1 Load data
1.1 Package
library(maptools)## Warning: package 'maptools' was built under R version 3.3.3## Loading required package: sp## Checking rgeos availability: TRUElibrary(maps)
library(mapdata)## Warning: package 'mapdata' was built under R version 3.3.3library(dismo)## Loading required package: rasterlibrary(raster)1.2 Load csv file
pachy <- read.csv("Herbario.csv")1.3 Download data 19 bioclimatic variables 2.5 minutes from worldclim
predStackModern <- getData("worldclim", var="bio", res=2.5, path="data/")
# set the coordinate reference system for WGS84
projection(predStackModern) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
# Crop database for region Neotropic
region = extent(-90,-30,-40,15) # define the extent neotropic
predStackModern_NEO = crop(predStackModern, region)
# Plot data climat
plot(predStackModern_NEO) 
2 Analysis
2.1 set output directory
dir.create('Output_modern_Maxent_NEO', recursive=TRUE) # creates new directory
dir.create('Output_modern_Maxent_WORLD', recursive=TRUE) 2.2 Train Maxent model with arguments
modelBasic_NEO <- maxent(
x=predStackModern_NEO,
p=pachy[ , c('Long', 'Lat')],
path=paste(getwd(), '/Output_modern_Maxent_NEO', sep=''),
args=c(
'randomtestpoints=30',
'betamultiplier=1',
'linear=true',
'quadratic=true',
'product=true',
'threshold=true',
'hinge=true',
'threads=2',
'responsecurves=true',
'jackknife=true',
'askoverwrite=false'
)
)2.3 look at model output (HTML page)
It will open HTML page in your browser (Firefox, Safari…)
modelBasic_NEO## class : MaxEnt
## variables: bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19Note that unlike running Maxent from the program the prediction raster is not automatically created, so we have to do that in the next step
3 Write prediction map
Note that you could save the prediction raster in any number of formats (see ?writeFormats), but GeoTiffs are small and can be read by ArcMap… ASCIIs can also be read by other programs but are large… the default raster format (GRD) sometimes can’t be read by ArcMap…if you don’t specifiy a file name then the results are written to R’s memory
mapModernBasic_NEO <- predict(
object=modelBasic_NEO,
x=predStackModern_NEO,
filename=paste(getwd(), '/Output_modern_Maxent_NEO/pachy_map_modern_neo', sep=''),
na.rm=TRUE,
format='GTiff',
overwrite=TRUE,
progress='text'
)##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
## Look at map
plot(mapModernBasic_NEO, main='Modern')
#Add species' records
points(pachy[ , c('Long', 'Lat')], cex=0.3)
4 FOR WORLDMAP
#MaxEnt
modelBasic <- maxent(
x=predStackModern,
p=pachy[ , c('Long', 'Lat')],
path=paste(getwd(), '/Output_modern_Maxent_WORLD', sep=''),
args=c(
'randomtestpoints=30',
'betamultiplier=1',
'linear=true',
'quadratic=true',
'product=true',
'threshold=true',
'hinge=true',
'threads=2',
'responsecurves=true',
'jackknife=true',
'askoverwrite=false'
)
)
# prediction
mapModernBasic <- predict(
object=modelBasic,
x=predStackModern,
filename=paste(getwd(), '/Output_modern_Maxent_WORLD/pachy_map_modern', sep=''),
na.rm=TRUE,
format='GTiff',
overwrite=TRUE,
progress='text'
)##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 6%
|
|===== | 8%
|
|====== | 9%
|
|======= | 11%
|
|======== | 12%
|
|========= | 14%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 18%
|
|============= | 20%
|
|============== | 21%
|
|=============== | 23%
|
|================ | 24%
|
|================= | 26%
|
|================== | 27%
|
|=================== | 29%
|
|==================== | 30%
|
|===================== | 32%
|
|====================== | 33%
|
|======================= | 35%
|
|======================== | 36%
|
|========================= | 38%
|
|========================== | 39%
|
|=========================== | 41%
|
|============================ | 42%
|
|============================= | 44%
|
|============================== | 45%
|
|=============================== | 47%
|
|================================ | 48%
|
|================================ | 50%
|
|================================= | 52%
|
|================================== | 53%
|
|=================================== | 55%
|
|==================================== | 56%
|
|===================================== | 58%
|
|====================================== | 59%
|
|======================================= | 61%
|
|======================================== | 62%
|
|========================================= | 64%
|
|========================================== | 65%
|
|=========================================== | 67%
|
|============================================ | 68%
|
|============================================= | 70%
|
|============================================== | 71%
|
|=============================================== | 73%
|
|================================================ | 74%
|
|================================================= | 76%
|
|================================================== | 77%
|
|=================================================== | 79%
|
|==================================================== | 80%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|=================================================================| 100%
## #Plot records
plot(mapModernBasic, main='Modern')
#Add species' records
points(pachy[ , c('Long', 'Lat')], cex=0.3)