Distribucion_ClaseConservacion
Jorge Lizarazo
3/7/2022
How to create a consensus Distribution model
It is to create a specie distribution model from different approaches. It means we are going to create some distribution analysis and then add those to only one.
Before starting: work space creation
Fist at all we must input and install the R packages needed, they are gonna be the next ones:
#install.packages(c("dismo","maptools", "sf", "raster","sp","dplyr", "kernlab", "rgdal", "randomForest", "kernlab"))note: Please to make an installing run once, no more times per computer!!!
Hence we are able to open or display our packages every time we want it
library(dismo)## Loading required package: raster## Loading required package: splibrary(maptools)## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.library(sp)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(raster)
library(rgdal)## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: D:/1. Mis Documentos/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: D:/1. Mis Documentos/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was D:/1. Mis Documentos/Documents/R/win-library/4.1/rgdal/projlibrary(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(randomForest)## randomForest 4.7-1## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combinelibrary(kernlab)##
## Attaching package: 'kernlab'## The following objects are masked from 'package:raster':
##
## buffer, rotatedSelect a work space or a folder for working on
You should type the folder address and take care about the order given with “/” icon.
setwd("F:/biato/wikiaves/Protocolo_distribucion_mapas/")
first_folder<-"F:/biato/wikiaves/Protocolo_distribucion_mapas/"note: there is an example. We would rather get you to have at list the last folder called as well the example, please. (“Protocolo_distribucion_mapas”)
Import data points from your data frame
You have should saved You have should saved your data frame as “.cvs” and with type CSV UTF (Comma delimited).
obs.data <- read.csv(file = "F:/biato/wikiaves/Protocolo_distribucion_mapas/Multuicolor.csv")note: there is an example “Multuicolor.csv”, but you data frame document should be called in a different way.
As you can see the data-frame from “GBif” has a lot of information that, although significant, are not relevant to our analysis. You can have a look from:
head(obs.data, 3)## gbifID datasetKey
## 1 3465991102 50c9509d-22c7-4a22-a47d-8c48425ef4a7
## 2 3457110124 50c9509d-22c7-4a22-a47d-8c48425ef4a7
## 3 3457106748 50c9509d-22c7-4a22-a47d-8c48425ef4a7
## occurrenceID kingdom phylum class
## 1 https://www.inaturalist.org/observations/105402005 Animalia Chordata Aves
## 2 https://www.inaturalist.org/observations/103054229 Animalia Chordata Aves
## 3 https://www.inaturalist.org/observations/103830088 Animalia Chordata Aves
## order family genus species
## 1 Passeriformes Thraupidae Chlorochrysa Chlorochrysa nitidissima
## 2 Passeriformes Thraupidae Chlorochrysa Chlorochrysa nitidissima
## 3 Passeriformes Thraupidae Chlorochrysa Chlorochrysa nitidissima
## infraspecificEpithet taxonRank scientificName
## 1 NA SPECIES Chlorochrysa nitidissima P.L.Sclater, 1874
## 2 NA SPECIES Chlorochrysa nitidissima P.L.Sclater, 1874
## 3 NA SPECIES Chlorochrysa nitidissima P.L.Sclater, 1874
## verbatimScientificName verbatimScientificNameAuthorship countryCode
## 1 Chlorochrysa nitidissima CO
## 2 Chlorochrysa nitidissima CO
## 3 Chlorochrysa nitidissima CO
## locality stateProvince occurrenceStatus individualCount
## 1 Valle del Cauca PRESENT NA
## 2 Caldas PRESENT NA
## 3 Valle del Cauca PRESENT NA
## publishingOrgKey decimalLatitude decimalLongitude
## 1 28eb1a3f-1c15-4a95-931a-4af90ecb574d 3.559347 -76.68557
## 2 28eb1a3f-1c15-4a95-931a-4af90ecb574d 5.372424 -75.96802
## 3 28eb1a3f-1c15-4a95-931a-4af90ecb574d 3.493145 -76.40396
## coordinateUncertaintyInMeters coordinatePrecision elevation elevationAccuracy
## 1 31422 NA NA NA
## 2 31384 NA NA NA
## 3 31422 NA NA NA
## depth depthAccuracy eventDate day month year taxonKey speciesKey
## 1 NA NA 2021-10-29T08:35:16 29 10 2021 5230597 5230597
## 2 NA NA 2021-10-17T13:55:00 17 10 2021 5230597 5230597
## 3 NA NA 2021-12-14T12:51:00 14 12 2021 5230597 5230597
## basisOfRecord institutionCode collectionCode catalogNumber recordNumber
## 1 HUMAN_OBSERVATION iNaturalist Observations 105402005
## 2 HUMAN_OBSERVATION iNaturalist Observations 103054229
## 3 HUMAN_OBSERVATION iNaturalist Observations 103830088
## identifiedBy dateIdentified license
## 1 Guillermo Nagy Aramacao Tours 2022-01-23T01:06:29 CC_BY_NC_4_0
## 2 David Monroy R 2021-12-14T17:57:49 CC_BY_NC_4_0
## 3 David J Barton 2021-12-28T17:28:26 CC_BY_NC_4_0
## rightsHolder recordedBy typeStatus
## 1 Guillermo Nagy Aramacao Tours Guillermo Nagy Aramacao Tours
## 2 David Monroy R David Monroy R
## 3 David J Barton David J Barton
## establishmentMeans lastInterpreted mediaType issue
## 1 2022-02-20T00:53:13.705Z COORDINATE_ROUNDED
## 2 2022-02-20T00:23:42.696Z StillImage COORDINATE_ROUNDED
## 3 2022-02-20T00:32:18.148Z StillImage COORDINATE_ROUNDEDnames(obs.data) ##number of colums - there are a lot of ## [1] "gbifID" "datasetKey"
## [3] "occurrenceID" "kingdom"
## [5] "phylum" "class"
## [7] "order" "family"
## [9] "genus" "species"
## [11] "infraspecificEpithet" "taxonRank"
## [13] "scientificName" "verbatimScientificName"
## [15] "verbatimScientificNameAuthorship" "countryCode"
## [17] "locality" "stateProvince"
## [19] "occurrenceStatus" "individualCount"
## [21] "publishingOrgKey" "decimalLatitude"
## [23] "decimalLongitude" "coordinateUncertaintyInMeters"
## [25] "coordinatePrecision" "elevation"
## [27] "elevationAccuracy" "depth"
## [29] "depthAccuracy" "eventDate"
## [31] "day" "month"
## [33] "year" "taxonKey"
## [35] "speciesKey" "basisOfRecord"
## [37] "institutionCode" "collectionCode"
## [39] "catalogNumber" "recordNumber"
## [41] "identifiedBy" "dateIdentified"
## [43] "license" "rightsHolder"
## [45] "recordedBy" "typeStatus"
## [47] "establishmentMeans" "lastInterpreted"
## [49] "mediaType" "issue"Clear and process our dataframe
We had better clear and delete some information for our data frame to avoid some mistakes and bad issues down to the code. We are going to take just the columns with geographic information owing to just one specie analysis. it is supposed to all coordinates been from one specie.
it is to select the columns 22 to 23, that are which has latitude and longitude respectively.
obs.data2 <- obs.data[,22:23]
head(obs.data2, 3)## decimalLatitude decimalLongitude
## 1 3.559347 -76.68557
## 2 5.372424 -75.96802
## 3 3.493145 -76.40396At the same time, we must delete the NAs data or empty boxes.
acgeo2 <- subset(obs.data2, !is.na(decimalLongitude) & !is.na(decimalLatitude))In order to adapt the dataframe for the function packages and simplify the code we will change the column order as well their names.
acgeo2<- select(acgeo2, decimalLongitude, decimalLatitude) ##reorder it
acgeo2<- rename(acgeo2, Long=decimalLongitude) #recall the column
acgeo2<- rename(acgeo2, Lat=decimalLatitude)#recall the columnnote: Long=decimalLongitude and Lat=decimalLatitud
Download the climatic layers from WWF data base and the politic composition raster of Colombia.
We are going to import some Bio Geographic rasters from internet to our exercise. There are some many ways to use more raster on our analysis and process but due to be a simple example we shall use just those.
files <- list.files(path=paste(system.file(package="dismo"),
'/ex', sep=''), pattern='grd', full.names=TRUE )
Bioma <- stack(files)note: Our object “Bioma” has all the imported raster which we are gonna use.
Downloading country political division
The “raster” package has some countries political division. In the Colombian case it has the country (layer 0), departments (layer 1), and regions (layer 2).
Colombia0 <- getData('GADM' , country="COL", level=0)
Colombia1 <- getData('GADM' , country="COL", level=1)
plot(Colombia0, main="Adm. Boundaries Colombia Level 0")
plot(Colombia1, main="Adm. Boundaries Colombia Level 1")
Keep proccesing it
Then we are gonna input our data points on the Colombian map downloaded
plot(Colombia0, xlim=c(-90,-60), ylim=c(-5,10), axes=TRUE, col="light yellow")
points(acgeo2$Lon, acgeo2$Lat, col='orange', pch=20, cex=0.75)
points(acgeo2$Lon, acgeo2$Lat, col='red', cex=0.75)
We need to cut the Bio climatic raster off by a specific geographic extension which must be given by us.
e <- extent(-81.841530, -66.87033, -4.228429 , 15.91247) #Colombia extent
Bioma <- crop(x = Bioma, y = e) ##this funtion cuts "Bioma" rasters by "e" extent
predictors2 <- stack(Bioma)
predictors2## class : RasterStack
## dimensions : 40, 30, 1200, 9 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -82, -67, -4, 16 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : bio1, bio12, bio16, bio17, bio5, bio6, bio7, bio8, biome
## min values : 77, 100, 80, 0, 135, 27, 64, 78, 1
## max values : 289, 7682, 2458, 1496, 361, 238, 157, 290, 14plot(predictors2)
Extract values for each point from our biogeographic rasters
We are going to extract values from our “Bioma” group of raster for each geographic point.
presvals <- extract(predictors2, acgeo2) #funtion to extract data
head(presvals, 3) #first registers ## bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## [1,] 224 2496 859 403 279 171 108 220 1
## [2,] 199 2312 725 409 256 146 110 197 1
## [3,] 183 1669 587 236 242 126 116 181 1tail(presvals, 3) #last registers## bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## [4664,] 202 2397 821 378 265 144 121 199 1
## [4665,] 192 2915 1014 281 247 130 117 190 1
## [4666,] 192 2915 1014 281 247 130 117 190 1#Take care if there are some 0 or NA values - pleas ask in that caseCreating adsence points
There are many ways to model a specie distribution, and yet most of them could be created by a blend of presence / absence information. To achieve our target we must create some pseudo-ramdom absence points.
set.seed(0)
backgr <- randomPoints(predictors2, 500)
absvals <- extract(predictors2, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
sdmdata[,'biome'] <- as.factor(sdmdata[,'biome'])
pred_nf <- dropLayer(predictors2, 'biome')
set.seed(0)
group <- kfold(acgeo2, 5)
pres_train <- acgeo2[group != 1, ]
pres_test <- acgeo2[group == 1, ]
########################
set.seed(100)
backg <- randomPoints(pred_nf, n=10000, ext=e, extf = 1.25)## Warning in randomPoints(pred_nf, n = 10000, ext = e, extf = 1.25): changed n to
## ncell of the mask (extent)## Warning in randomPoints(pred_nf, n = 10000, ext = e, extf = 1.25): generated
## random points = 0.654166666666667 times requested numbercolnames(backg) = c('Lon', 'Lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
r <- raster(pred_nf, 1)
plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE)
plot(e, add=TRUE, col='red', lwd=2)
points(backg_train, pch='-', cex=0.5, col='yellow')
points(backg_test, pch='-', cex=0.5, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')
Modeling specie distribution
Now we have all the stuff to make some species distribution models. We are going to set one by one and at the end make an unified model.
##Bioclim classic method (or climate-envelope-model)
bc <- bioclim(pred_nf, pres_train)
plot(bc, a=1, b=2, p=0.85)
ef <- evaluate(pres_test, backg_test, bc, pred_nf)
ef## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9854759
## cor : 0.5315259
## max TPR+TNR at : 0.0098116tr <- threshold(ef, 'spec_sens')
tr## [1] 0.0098116pb <- predict(pred_nf, bc, ext=e, progress='')
pb## class : RasterLayer
## dimensions : 40, 30, 1200 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -82, -67, -4, 16 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 0.5119207 (min, max)###### Grafica
par(mfrow=c(1,2))
plot(pb, main='Bioclim, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
plot(pb > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
Domain method
dm <- domain(pred_nf, pres_train)
ef <- evaluate(pres_test, backg_test, dm, pred_nf)
ef## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9456278
## cor : 0.7024954
## max TPR+TNR at : 0.7821892##Graphic
pd = predict(pred_nf, dm, ext=e, progress='')
par(mfrow=c(1,2))
plot(pd, main='Domain, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(ef, 'spec_sens')
plot(pd > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='*')
Mahalanobis distance
mm <- mahal(pred_nf, pres_train)
ef <- evaluate(pres_test, backg_test, mm, pred_nf)
ef## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9680914
## cor : 0.3650916
## max TPR+TNR at : 0.9999pm = predict(pred_nf, mm, ext=e, progress='')
par(mfrow=c(1,2))
pm[pm < -10] <- -10
##Graphic
plot(pm, main='Mahalanobis distance')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(ef, 'spec_sens')
plot(pm > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
Classical regression models
data3 <- backg_train
colnames(data3) <- colnames(pres_train)
rbind(pres_train, data3)
backg_train<-data3
train <- rbind(pres_train, backg_train)
pb_train <- c(rep(1, nrow(pres_train)), rep(0, nrow(backg_train)))
envtrain <- extract(predictors2, train)
envtrain <- data.frame( cbind(pa=pb_train, envtrain) )
envtrain[,'biome'] = factor(envtrain[,'biome'], levels=1:14)head(envtrain)## pa bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## 1 1 224 2496 859 403 279 171 108 220 1
## 2 1 199 2312 725 409 256 146 110 197 1
## 3 1 183 1669 587 236 242 126 116 181 1
## 4 1 224 2496 859 403 279 171 108 220 1
## 5 1 224 2496 859 403 279 171 108 220 1
## 6 1 202 2397 821 378 265 144 121 199 1testpres <- data.frame( extract(predictors2, pres_test) )
testbackg <- data.frame( extract(predictors2, backg_test) )
testpres[ ,'biome'] = factor(testpres[ ,'biome'], levels=1:14)
testbackg[ ,'biome'] = factor(testbackg[ ,'biome'], levels=1:14)Generalized Linear Models
gm1 <- glm(pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
family = binomial(link = "logit"), data=envtrain)
coef(gm1)## (Intercept) bio1 bio5 bio6 bio7 bio8
## 23.860298281 1.800882008 -0.552967161 -0.215889158 0.126127559 -1.125627220
## bio12 bio16 bio17
## -0.002742331 0.004338260 0.010615441gm2 <- glm(pa ~ bio1+bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
family = gaussian(link = "identity"), data=envtrain)
evaluate(testpres, testbackg, gm1)## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.8808139
## cor : 0.6409888
## max TPR+TNR at : 2.086457ge2 <- evaluate(testpres, testbackg, gm2)
ge2## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.8892041
## cor : 0.677396
## max TPR+TNR at : 0.8190103pg <- predict(predictors2, gm2, ext=e)
par(mfrow=c(1,2))
plot(pg, main='GLM/gaussian, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(ge2, 'spec_sens')
plot(pg > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)
Machine learning methods
to archive this step you must install a new package that is called “rJava”.
#install.packages("rJava") #to run this step only one time per computer
library(rJava)
###just in case you have an error please run the next few codes to update your "Java" full plugins. at the same time if the issue going on creating an error you should download a new "Java" program from your web browser.Then you are able to continue installing and libraring your "rJava" package.
#Sys.setenv(JAVA_HOME='C:\\Program Files\\Java\\jre8') # for 64-bit version
#Sys.setenv(JAVA_HOME='C:\\Program Files (x86)\\Java\\jre8') # for 32-bit versionHence we can go forward…
xm <- maxent(predictors2, pres_train, factors='biome')## Warning in .local(x, p, ...): only got:756random background point values; Small
## exent? Or is there a layer with many NA values?## This is MaxEnt version 3.4.3plot(xm)
response(xm)
ef <- evaluate(pres_test, backg_test, xm, predictors2)
ef## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9388624
## cor : 0.7255079
## max TPR+TNR at : 0.5400576px <- predict(predictors2, xm, ext=e, progress='')
par(mfrow=c(1,2))
## Grafica
plot(px, main='Maxent, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(ef, 'spec_sens')
plot(px > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
Random Forest
It is a good approach, if your specie actually is a restricted one.
model <- pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17
rf1 <- randomForest(model, data=envtrain)## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?model <- factor(pa) ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17
rf2 <- randomForest(model, data=envtrain)
rf3 <- randomForest(envtrain[,1:8], factor(pb_train))
erf <- evaluate(testpres, testbackg, rf1)
erf## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9749899
## cor : 0.955147
## max TPR+TNR at : 0.9220572pr <- predict(predictors2, rf1, ext=e)
##Graphic
par(mfrow=c(1,2))
plot(pr, main='Random Forest, regression')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(erf, 'spec_sens')
plot(pr > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)
Support Vector Machines method
svm <- ksvm(pa ~ bio1+bio5+bio6+bio7+bio8+bio12+bio16+bio17, data=envtrain)
esv <- evaluate(testpres, testbackg, svm)
esv## class : ModelEvaluation
## n presences : 933
## n absences : 157
## AUC : 0.9683816
## cor : 0.9379215
## max TPR+TNR at : 0.9432463ps <- predict(predictors2, svm, ext=e)
### Graphic
par(mfrow=c(1,2))
plot(ps, main='Support Vector Machine')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(esv, 'spec_sens')
plot(ps > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)
Creating a consensus method
Taking our models by different methods we can create a consensus model
models <- stack(pb, pd, pm, pg, pr, ps)
names(models) <- c("bioclim", "domain", "mahal", "glm", "rf", "svm")
##if your specie dose not inhabit forest please take out "rf" method from the last step.
plot(models)
m <- mean(models)
m<- mask(m,Colombia0 )
# Graphic
plot(m, main='average score')
points(pres_train, pch='.')
points(backg_train, pch='-', cex=0.25)
plot(Colombia0, add=TRUE, border='dark grey')
On the other hand, if you specie is restricted one, please use forward the next code:
auc <- sapply(list(ge2, erf, esv), function(x) x@auc)
w <- (auc-0.5)^2
m2 <- weighted.mean( models[[c("glm", "rf", "svm")]], w)
plot(m2, main='weighted mean of three models')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)
note: You must take care and use the object “m” or “m2” on the next field depends on your specie type.
Exporting specie distribution model given and saved as a polygon
First things first, at that point We would like to save the information in a new folder to process the next steps in our procedure.We could called “distribucion” folder.
file <- "distribucion"
dir.create(file.path(first_folder, file )) ## Warning in dir.create(file.path(first_folder, file)): 'F:
## \biato\wikiaves\Protocolo_distribucion_mapas\\distribucion' already exists ##Create a new folder called "distribucion"
##You shall see it on your folders
setwd("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion")Now, there are many types of Geographic documents, but in our case the best way is a Polygon raster, and yet it must be the specie presence above 80% probability.
setwd("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion")
pol <- rasterToPolygons(m>0.8,function(x) x == 1,dissolve=T) ##depends on your specie you would increase or decrease the porcentage ## Loading required namespace: rgeos#plot(Colombia0, axes=TRUE,col="light yellow", bg="light blue")
#plot(pol, add=T, col="green") ## it is just in case you wanna look this out
raster::shapefile(pol,
"F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion/multicolor_fin.shp") ### it import the final polygon raster in created folderPlotting our final specie distribution to look this out in detail
We have rather take a look and plot it in a better way, then we could have a dynamic view of it.
It is going to be from a package that get us a dynamic visualization. it is called “leaflet” and over there you can move, scroll, and zoom your polygon distribution over a mapping as good as Google map. So, we need to change the work space and install some new packages.
#install.packages("htmlwidgets", "webshot", "leaflet", "mapview")
#the installation funtion must be runed just ones per computer, please
library(htmlwidgets)
library(webshot)
library(leaflet)
library(mapview)
getwd()The next step is to create a new function which would be able to extract the average coordinates from our polygon raster. It would be useful for plotting the distribution over the accurate coordinates and zooming.
library(sf)
wqp_sites <- st_read("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion")## Reading layer `multicolor_fin' from data source
## `F:\biato\wikiaves\Protocolo_distribucion_mapas\distribucion'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77 ymin: 3.5 xmax: -75 ymax: 7
## Geodetic CRS: GCS_unknownnhd_wms_url <- "https://server.arcgisonline.com/ArcGIS/rest/services/World_Physical_Map/MapServer/tile/{z}/{y}/{x}"
setwd("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion")
extractCoords <- function(sp.df)
{
results <- list()
for(i in 1:length(sp.df@polygons[[1]]@Polygons))
{
results[[i]] <- sp.df@polygons[[1]]@Polygons[[i]]@coords
}
results <- Reduce(rbind, results)
results
}Hence, we will able to keep plotting our distribution
setwd("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion")
nhd_wms_url <- "https://server.arcgisonline.com/ArcGIS/rest/services/World_Physical_Map/MapServer/tile/{z}/{y}/{x}"
for(i in list.files()[grep(list.files(), pattern="[.]shp$")]){
dc=readOGR(i)
ext=extractCoords(dc) #Estraigo las coordenadas del shape file. con una fución de arriba
df<-data.frame(longmax=max(ext[,1]),
longmin=min(ext[,1]), #Saco el minimo y maximo de long y lat para
latmax=max(ext[,2]), #que después le pueda decir que la imagen quede con zoom
latmin=min(ext[,2])) #y centrada donde hay coordenadas. todo automaticamente
plotee<-print(leaflet() %>%
addTiles() %>%
fitBounds(df$longmin,df$latmin,df$longmax,df$latmax) %>%
addTiles() %>%
addMapPane("background_map", zIndex = 410) %>% # Level 1: bottom
addMapPane("polygons", zIndex = 420) %>% # Level 2: middle
addMapPane("labels", zIndex = 430) %>%
addWMSTiles(nhd_wms_url, layers = "2",
options = pathOptions(pane = "background_map") ) %>%
addPolygons(data = dc,
color = "71E775",
fillColor = "#71E775",
fillOpacity = .7,
weight = 1,
options = pathOptions(pane = "polygons"))%>%
addProviderTiles("Stamen.TonerHybrid",
options = pathOptions(pane = "labels"))) # este es pa que salgan las ciudades
mapshot(plotee, file = paste0("F:/biato/wikiaves/Protocolo_distribucion_mapas/distribucion",sub(".shp",".png",i)),#Este es para
vwidth = 1100, vheight = 1300) #guardar las imagenes.
}OGR data source with driver: ESRI Shapefile Source: “F:_distribucion_mapas_fin.shp”, layer: “multicolor_fin” with 1 features It has 1 fields
plotee