# necessary packages
library(sp)     # provides classes for spatial data in R
library(raster) # provides classes and methods for raster datasets

# necessary - but attention for Mac users!
library(rgdal)  # interface to the Geospatial Data Abstraction Library to read and write different spatial file formats

# usefull packages
library(maptools)   # some specific spatial methods and conversion methods for spatial data 
library(data.table)
library(fields)     # contains methods for spatial interpolation and statistics but also nice color palettes fro mapping
library(RColorBrewer) # color palettes for mapping
library(marxan)
library(rasterVis)
library(viridis) # color palette
library(RColorBrewer) # color palette
#library(gstat)

library(stringr)
library(dplyr)
# library(ggplot2)
library(gridExtra)
library (readxl)
library (sf)
library(tmap) # mapas
library(velox)
library(sp)
library(corrplot) # plot correlations
library(kableExtra) # tables

Introducción

Servicios ecosistemicos y especies Selección de áreas prioritarias con Marxan. Teniendo en cuenta el la distribucion de especies amenazadas (73 sp con mapas en biomodelos) y Mamiferos. Tambien se tuvo en cuenta la red actual de areas protegidas y un costo representado por la aptitud agricola de UPRA.

Ver Wilson et al 2011.
Ver Cuesta et al 2017

Método breve

Servicios: modelados con Invest Mauro

Carbono

Nutrientes

Agua

## Set coodenadas
crs.geo <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")    # geographical, datum WGS84

####################
#### Read Data
####################

path_offi <- "D:/BoxFiles/"
# path_offi <- "C:/Users/diego.lizcano/"
  
### limite colombia
col_limit1 <- readOGR(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/shp",sep=""),
                layer="mask_co", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp", layer: "mask_co"
## with 142 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID FID_siluet ID GRIDCODE FID_tapar Id_1

col_limit <- spTransform(col_limit1, crs.geo) # transform to geografico


#### read Carbon
carbon1 <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/Carbon/carbon_2/tot_c_cur_1k.tif", sep=""))

#### read Nutrientes
nutrients1 <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/nutrients/n_export_1k.tif", sep="")) 
### Aporte Nitrogeno
nutrients2 <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/nutrients/retencion_nutrientes/nutrient_ret.tif", sep="")) # retencion nut

sediments1 <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/sediments/sediment_ret1.tif", sep=""))

agua1 <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/water_yield/wyield_vol.tif", sep=""))

nutrients1 <- projectRaster(nutrients2, crs=crs.geo) # transform to geografico
sediments <- projectRaster(sediments1, nutrients1) # transform to geografico
carbon <-  projectRaster(carbon1, nutrients1) # reproject to fix 2 pixel desfase

water <- projectRaster(agua1, nutrients1) # reproject to fix 2 pixel desfase

# carbon <- projectRaster(carbon, crs=crs.geo)
   
servicios1 <- stack(carbon, nutrients1, sediments, water) 
names(servicios1) <- c("Carbon", "Nutrients", "Sediments", "Water")



servicios <- projectRaster(servicios1, crs=crs.geo)
plot(servicios)

# set the window extent 
w <- as(extent(servicios), "SpatialPolygons") # make extent a polygon
proj4string(w) <- crs.geo # put georef
w$id <- 1 # put value

w <- rasterize(w, carbon, 'id', fun=min) # get the bounding box of carbon

Selección de especies

Se apilaron los modelos de distribución de Biomodelos del IAvH, seleccionando los grupos: Especies amenazadas, Mamiferos, aves reptiles, y anfibios. Los modelos apilados corelacionaron espacialmente con los servicios ecosistemicos

Modelos de distribución

Se descargaron 73 especies amenazadas y 116 mamiferos xx aves, xx reptiles xx anfibios

Load vertebrate layers

#### 
mammal <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/biodiversity/Mammals_no_crop.tif", sep=""))

#### 
reptiles <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/biodiversity/Reptilia_no_crop.tif", sep="")) 
### 
aves <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/biodiversity/Aves_no_crop.tif", sep="")) # retencion nut

anfibios <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/biodiversity/Amphibia_no_crop.tif", sep=""))

amenazadas <- raster(paste(path_offi, "Box Sync/CodigoR/Servicios_Biodiversidad/raster/biodiversity/Amenazadas_no_crop.tif", sep=""))


# #re-define the projection and extent using nutrients1
mammal <- projectRaster(mammal, nutrients1) # transform to geografico
reptiles <- projectRaster(reptiles, nutrients1) # transform to geografico
aves <-  projectRaster(aves, nutrients1)
anfibios <-  projectRaster(anfibios, nutrients1)
amenazadas <-  projectRaster(amenazadas, nutrients1)

# Crea un stack
especies <- stack(mammal, reptiles, aves, anfibios, amenazadas)
names(especies) <- c("Mammal", "Reptile", "Bird", "Amphibia", "Threatened")

plot(especies)
plot(col_limit, add=T)

my_colors <- viridis(9)

#### see as thematic map using tmap package
tm_shape(especies) + 
    tm_raster(c("Mammal", "Reptile", "Bird", "Amphibia", "Threatened"), 
              breaks=c(-Inf, 1, 10, 20, 50, 100, 200, 300, 400, 600),  
        palette = my_colors,  auto.palette.mapping=FALSE) +  
  tm_shape(col_limit) + 
  tm_borders() +
  tm_legend(text.size=0.8,
    # title.size=1.2,
    position = c("left", "bottom"),
    panel.labels=c("Mammal", "Reptile", "Bird", "Amphibia", "Threatened"),
    bg.color = "white", 
    bg.alpha=.25, 
    frame="gray50", 
    height=.6, 
    hist.width=.2,
    hist.height=.2, 
    hist.bg.color="gray60", 
    hist.bg.alpha=.5)

Results

Correlations

# load ecoregions
Amazonia<-readOGR(paste(path_offi,"Box Sync/CodigoR/Servicios_Biodiversidad/shp/Amazonia", sep=""),
                layer="408_Amazonia", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\Amazonia", layer: "408_Amazonia"
## with 1 features
## It has 9 fields
## Integer64 fields read as strings:  COUNT Z_OF_SYS

Caribbean<-readOGR(paste(path_offi,"Box Sync/CodigoR/Servicios_Biodiversidad/shp/Caribbean", sep=""),
                layer="411_Caribbean", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\Caribbean", layer: "411_Caribbean"
## with 1 features
## It has 9 fields
## Integer64 fields read as strings:  COUNT Z_OF_SYS

Meso_American<-readOGR(paste(path_offi,"Box Sync/CodigoR/Servicios_Biodiversidad/shp/Meso_american_lowland", sep=""),
                layer="420_Meso-American_Lowland_Forest", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\Meso_american_lowland", layer: "420_Meso-American_Lowland_Forest"
## with 1 features
## It has 9 fields
## Integer64 fields read as strings:  COUNT Z_OF_SYS

Andes<-readOGR(paste(path_offi,"Box Sync/CodigoR/Servicios_Biodiversidad/shp/North_central_moist_andes", sep=""),
                layer="409_North-Central_Moist_Andes", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\North_central_moist_andes", layer: "409_North-Central_Moist_Andes"
## with 1 features
## It has 9 fields
## Integer64 fields read as strings:  COUNT Z_OF_SYS

Orinoquia<-readOGR(paste(path_offi,"Box Sync/CodigoR/Servicios_Biodiversidad/shp/Orinoquia", sep=""),
                layer="405_Orinoquia", stringsAsFactors=FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\Orinoquia", layer: "405_Orinoquia"
## with 1 features
## It has 9 fields
## Integer64 fields read as strings:  COUNT Z_OF_SYS

##### Corrrelations
#### Function to get probability
cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], na.omit=T ,...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
    colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
    p.mat
} # end of function



# crea stack con todo
all_stack <- addLayer(servicios1, especies)  

# Mask by amazonas
all_stack_masked_a <- mask(all_stack, Amazonia)
# matrix of the p-value of the correlation
p.mat_a <- cor.mtest(all_stack_masked_a)

jnk=layerStats(all_stack_masked_a, 'pearson', na.rm=T)
corr_matrix_a=jnk$'pearson correlation coefficient'
par(mfrow=c(1,2)) 
plot(col_limit)
plot(Amazonia, col="grey", add=T) 
title("Amazonia", line = -1)
corrplot(corr_matrix_a, type = "upper", method = "circle", 
         p.mat = p.mat_a, sig.level = 0.01)

# get tables
kable(corr_matrix_a, format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	Carbon	Nutrients	Sediments	Water	Mammal	Reptile	Bird	Amphibia	Threatened
Carbon	1.0000000	0.6959749	-0.3457926	0.1691980	-0.1820663	-0.0782093	-0.2136615	-0.1715534	-0.1342105
Nutrients	0.6959749	1.0000000	-0.3874752	0.1601722	-0.1742258	-0.0679759	-0.2056112	-0.1561962	-0.1264445
Sediments	-0.3457926	-0.3874752	1.0000000	-0.0457698	0.2100047	0.0321206	0.2196295	0.1385425	0.1532857
Water	0.1691980	0.1601722	-0.0457698	1.0000000	0.0029241	-0.1174603	-0.0729246	-0.0743108	-0.0506639
Mammal	-0.1820663	-0.1742258	0.2100047	0.0029241	1.0000000	0.7956714	0.8232246	0.9317667	0.8073344
Reptile	-0.0782093	-0.0679759	0.0321206	-0.1174603	0.7956714	1.0000000	0.8178706	0.8509066	0.8374619
Bird	-0.2136615	-0.2056112	0.2196295	-0.0729246	0.8232246	0.8178706	1.0000000	0.8162272	0.9168312
Amphibia	-0.1715534	-0.1561962	0.1385425	-0.0743108	0.9317667	0.8509066	0.8162272	1.0000000	0.7707557
Threatened	-0.1342105	-0.1264445	0.1532857	-0.0506639	0.8073344	0.8374619	0.9168312	0.7707557	1.0000000

# Mask by Caribean
all_stack_masked_c <- mask(all_stack, Caribbean)
# matrix of the p-value of the correlation
p.mat_c <- cor.mtest(all_stack_masked_c)

jnk=layerStats(all_stack_masked_c, 'pearson', na.rm=T)
corr_matrix_c=jnk$'pearson correlation coefficient'
par(mfrow=c(1,2)) 
plot(col_limit)
plot(Caribbean,  col="grey", add=T )
title("Caribbean", line = -1)
corrplot(corr_matrix_c, type = "upper", method = "circle",
         p.mat = p.mat_c, sig.level = 0.01)

kable(corr_matrix_c, format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	Carbon	Nutrients	Sediments	Water	Mammal	Reptile	Bird	Amphibia	Threatened
Carbon	1.0000000	0.4486914	0.2459836	0.2209631	0.0516946	-0.2203186	0.2258081	-0.1405714	0.3494772
Nutrients	0.4486914	1.0000000	0.2355710	0.2704853	0.1215836	-0.1307965	0.2499940	-0.0350756	0.3119661
Sediments	0.2459836	0.2355710	1.0000000	0.6818269	0.1935819	-0.1896959	0.1817719	-0.0307003	0.5387033
Water	0.2209631	0.2704853	0.6818269	1.0000000	0.5231647	0.1798495	0.3902003	0.3875528	0.6104257
Mammal	0.0516946	0.1215836	0.1935819	0.5231647	1.0000000	0.6434260	0.3736993	0.6467665	0.3156902
Reptile	-0.2203186	-0.1307965	-0.1896959	0.1798495	0.6434260	1.0000000	-0.1031572	0.8668685	-0.3082301
Bird	0.2258081	0.2499940	0.1817719	0.3902003	0.3736993	-0.1031572	1.0000000	0.1873139	0.6157180
Amphibia	-0.1405714	-0.0350756	-0.0307003	0.3875528	0.6467665	0.8668685	0.1873139	1.0000000	-0.1127891
Threatened	0.3494772	0.3119661	0.5387033	0.6104257	0.3156902	-0.3082301	0.6157180	-0.1127891	1.0000000

# Mask by Meso_American
all_stack_masked_m <- mask(all_stack, Meso_American)
# matrix of the p-value of the correlation
p.mat_m <- cor.mtest(all_stack_masked_m)

jnk=layerStats(all_stack_masked_m, 'pearson', na.rm=T)
corr_matrix_m=jnk$'pearson correlation coefficient'
par(mfrow=c(1,2)) 
plot(col_limit)
plot(Meso_American, col="grey", add=T )
title("Meso_American", line = -1)
corrplot(corr_matrix_m, type = "upper", method = "circle", 
         p.mat = p.mat_m, sig.level = 0.01)

kable(corr_matrix_m, format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	Carbon	Nutrients	Sediments	Water	Mammal	Reptile	Bird	Amphibia	Threatened
Carbon	1.0000000	0.7457283	0.1959647	0.5722948	-0.0433094	-0.2073358	-0.2034249	-0.1612182	0.1143046
Nutrients	0.7457283	1.0000000	0.1614918	0.5182750	-0.0324159	-0.1658559	-0.1869700	-0.1211938	0.1028128
Sediments	0.1959647	0.1614918	1.0000000	0.1709904	0.1733565	0.0263272	-0.0817930	0.0960048	0.1594391
Water	0.5722948	0.5182750	0.1709904	1.0000000	-0.1511210	-0.3536843	-0.5068695	-0.2548885	0.0124271
Mammal	-0.0433094	-0.0324159	0.1733565	-0.1511210	1.0000000	0.8076963	0.6411362	0.8440115	0.8221579
Reptile	-0.2073358	-0.1658559	0.0263272	-0.3536843	0.8076963	1.0000000	0.5698629	0.9244844	0.5708307
Bird	-0.2034249	-0.1869700	-0.0817930	-0.5068695	0.6411362	0.5698629	1.0000000	0.5296217	0.7102935
Amphibia	-0.1612182	-0.1211938	0.0960048	-0.2548885	0.8440115	0.9244844	0.5296217	1.0000000	0.6097827
Threatened	0.1143046	0.1028128	0.1594391	0.0124271	0.8221579	0.5708307	0.7102935	0.6097827	1.0000000

# Mask by Andes
all_stack_masked_an <- mask(all_stack, Andes)
# matrix of the p-value of the correlation
p.mat_an <- cor.mtest(all_stack_masked_an)

jnk=layerStats(all_stack_masked_an, 'pearson', na.rm=T)
corr_matrix_an=jnk$'pearson correlation coefficient'
par(mfrow=c(1,2)) 
plot(col_limit)
plot(Andes, col="grey", add=T )
title("Andes", line = -1)
corrplot(corr_matrix_an, type = "upper", method = "circle",
         p.mat = p.mat_an, sig.level = 0.01)

kable(corr_matrix_an, format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	Carbon	Nutrients	Sediments	Water	Mammal	Reptile	Bird	Amphibia	Threatened
Carbon	1.0000000	0.7684310	0.2386746	0.3760061	0.0571707	-0.0421356	-0.0018815	0.0071693	0.1593514
Nutrients	0.7684310	1.0000000	0.2564793	0.3798009	0.0091817	-0.0803514	-0.0397636	-0.0226970	0.1591400
Sediments	0.2386746	0.2564793	1.0000000	0.5479447	0.1781147	0.0542397	0.0480534	0.1343344	0.2022619
Water	0.3760061	0.3798009	0.5479447	1.0000000	0.3403386	0.2319615	0.0018688	0.3338561	0.2451139
Mammal	0.0571707	0.0091817	0.1781147	0.3403386	1.0000000	0.8445903	0.6187383	0.8587301	0.6794045
Reptile	-0.0421356	-0.0803514	0.0542397	0.2319615	0.8445903	1.0000000	0.3699837	0.8975206	0.3357763
Bird	-0.0018815	-0.0397636	0.0480534	0.0018688	0.6187383	0.3699837	1.0000000	0.4606414	0.7391988
Amphibia	0.0071693	-0.0226970	0.1343344	0.3338561	0.8587301	0.8975206	0.4606414	1.0000000	0.5505224
Threatened	0.1593514	0.1591400	0.2022619	0.2451139	0.6794045	0.3357763	0.7391988	0.5505224	1.0000000

# Mask by Orinoquia
all_stack_masked_o <- mask(all_stack, Orinoquia)
# matrix of the p-value of the correlation
p.mat_o <- cor.mtest(all_stack_masked_o)

jnk=layerStats(all_stack_masked_o, 'pearson', na.rm=T)
corr_matrix_o=jnk$'pearson correlation coefficient'
par(mfrow=c(1,2)) 
plot(col_limit)
plot(Orinoquia, col="grey", add=T )
title("Orinoquia", line = -1)
corrplot(corr_matrix_o, type = "upper", method = "circle",
         p.mat = p.mat_o, sig.level = 0.01)

# ecoregions$NS_DIVISIO has names
# 
# names(sp_amen_richnes) <- "Species"    
# ### correlations     
# all_stack <- addLayer(servicios1, sp_amen_richnes)   
# all_stack <- addLayer(all_stack, mammal_stack) 
# library(corrplot)
# jnk=layerStats(all_stack, 'pearson', na.rm=T)
# corr_matrix=jnk$'pearson correlation coefficient'
# corrplot(corr_matrix, type = "upper", method = "circle")

Statistics of correlations

kable(corr_matrix_o, format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

	Carbon	Nutrients	Sediments	Water	Mammal	Reptile	Bird	Amphibia	Threatened
Carbon	1.0000000	0.3370423	0.0286969	0.0502950	0.1049423	0.0799174	0.0494925	0.0426281	0.1209600
Nutrients	0.3370423	1.0000000	-0.1232427	-0.0809510	-0.0973369	-0.1294697	-0.0819302	-0.1700609	-0.0496550
Sediments	0.0286969	-0.1232427	1.0000000	0.1163070	0.4975159	0.4037481	0.4668249	0.5203646	0.4210377
Water	0.0502950	-0.0809510	0.1163070	1.0000000	0.4354258	0.4825927	0.1757094	0.4086294	0.5786776
Mammal	0.1049423	-0.0973369	0.4975159	0.4354258	1.0000000	0.8316978	0.7451743	0.9193044	0.8432187
Reptile	0.0799174	-0.1294697	0.4037481	0.4825927	0.8316978	1.0000000	0.5515067	0.8238283	0.7463610
Bird	0.0494925	-0.0819302	0.4668249	0.1757094	0.7451743	0.5515067	1.0000000	0.7353320	0.6826957
Amphibia	0.0426281	-0.1700609	0.5203646	0.4086294	0.9193044	0.8238283	0.7353320	1.0000000	0.7696441
Threatened	0.1209600	-0.0496550	0.4210377	0.5786776	0.8432187	0.7463610	0.6826957	0.7696441	1.0000000

Map correlations

# 
# RasterStack or Brick objects
  Amaz_mamif_sedimen <- corLocal(all_stack_masked_a[[5]],
                 all_stack_masked_a[[3]],
                ngb=11, method=c("pearson"))
  
  Amaz_ave_sedimen <-  corLocal(all_stack_masked_a[[7]],
                 all_stack_masked_a[[3]],
                ngb=11, method=c("pearson"))
  
  Amaz_amenaz_sedimen <-  corLocal(all_stack_masked_a[[9]],
                 all_stack_masked_a[[3]],
                ngb=11, method=c("pearson"))
  
  Caribe_agua_amenaz <-  corLocal(all_stack_masked_c[[4]],
                 all_stack_masked_c[[9]],
                ngb=11, method=c("pearson"))
    
  Caribe_sedimen_amenaz <-  corLocal(all_stack_masked_c[[3]],
                 all_stack_masked_c[[9]],
                ngb=11, method=c("pearson")) 
  
  Caribe_agua_reptil <-  corLocal(all_stack_masked_c[[4]],
                 all_stack_masked_c[[6]],
                ngb=11, method=c("pearson")) 
    
  Meso_sediment_amenaz <-  corLocal(all_stack_masked_m[[3]],
                 all_stack_masked_m[[9]],
                ngb=11, method=c("pearson")) 
    
  Meso_carbon_amenaza <-  corLocal(all_stack_masked_m[[1]],
                 all_stack_masked_m[[9]],
                ngb=11, method=c("pearson")) 
    
  Andes_agua_amenaz <- corLocal(all_stack_masked_an[[4]],
                 all_stack_masked_an[[9]],
                ngb=11, method=c("pearson"))  
  
  Andes_agua_anfi <- corLocal(all_stack_masked_an[[4]],
                 all_stack_masked_an[[8]],
                ngb=11, method=c("pearson"))  
  
  Orinoq_agua_anfi <- corLocal(all_stack_masked_o[[4]],
                 all_stack_masked_o[[9]],
                ngb=11, method=c("pearson")) 
  
  Orinoq_Sediments_anfi <- corLocal(all_stack_masked_o[[3]],
                 all_stack_masked_o[[9]],
                ngb=11, method=c("pearson")) 
    
    
 
  
# stack  
corr_stack <- stack(Amaz_mamif_sedimen, 
                    Amaz_ave_sedimen,
                    Caribe_agua_amenaz,
                    Caribe_sedimen_amenaz,
                    Meso_sediment_amenaz,
                    Meso_carbon_amenaza,
                    Andes_agua_amenaz,
                    Andes_agua_anfi, 
                    Orinoq_agua_anfi,
                    Orinoq_Sediments_anfi)
# names
names(corr_stack)<-c("Mamiferos_Sedimento", 
                     "Aves_Sedimento",
                     "Amenazadas_Agua1", 
                     "Amenazadas_Sedimento1", 
                     "Amenazadas_Sedimento2",
                     "Amenazadas_Carbono", 
                     "Amenazadas_Agua2",
                     "Anfibios_Agua1",
                     "Anfibios_Agua2",
                     "Anfibios_Sedimento")
# plot(corr_stack, col=brewer.pal('RdBu', n=20))
my_color2 <- brewer.pal('RdBu', n=10)

#### see as thematic map using tmap package
tm_shape(corr_stack) + 
    tm_raster(c("Mamiferos_Sedimento", 
                     "Aves_Sedimento",
                     "Amenazadas_Agua1", 
                     "Amenazadas_Sedimento1", 
                     "Amenazadas_Sedimento2",
                     "Amenazadas_Carbono", 
                     "Amenazadas_Agua2",
                     "Anfibios_Agua1",
                     "Anfibios_Agua2",
                     "Anfibios_Sedimento"), 
              breaks=c(-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1),  
        palette = my_color2) +  
        # auto.palette.mapping=FALSE) +  
  tm_shape(col_limit) + 
  tm_borders() +
  # tm_legend(text.size=0.5,
  #   title.size=0.7,
  #   position = c("left", "bottom"),
  #   # panel.labels=c("Mamiferos_Sedimento", 
  #   #                  "Aves_Sedimento",
  #   #                  "Amenazadas_Agua1", 
  #   #                  "Amenazadas_Sedimento1", 
  #   #                  "Amenazadas_Sedimento2",
  #   #                  "Amenazadas_Carbono", 
  #   #                  "Amenazadas_Agua2",
  #   #                  "Anfibios_Agua1",
  #   #                  "Anfibios_Agua2",
  #   #                  "Anfibios_Sedimento"),
  #   bg.color = "white", 
  #   bg.alpha=.25, 
  #   frame="gray50", 
  #   height=.6, 
  #   hist.width=.2,
  #   hist.height=.2, 
  #   hist.bg.color="gray60", 
  #   hist.bg.alpha=.5)

  tm_layout(
     panel.labels=c("Mamiferos_Sedimento",
                      "Aves_Sedimento",
                     "Amenazadas_Agua1",
                      "Amenazadas_Sedimento1",
                      "Amenazadas_Sedimento2",
                      "Amenazadas_Carbono",
                      "Amenazadas_Agua2",
                     "Anfibios_Agua1",
                     "Anfibios_Agua2",
                      "Anfibios_Sedimento"),
    legend.show = FALSE)

p <- levelplot(corr_stack, par.settings = RdBuTheme)
p + layer(sp.lines(col_limit, lwd=0.8, col='darkgray'))

# 
# rm(corr_stack, yc_r, yn_r, ys_r, ys_m )

Congruencia Espacial Entre Biodiversidad y Servicios Ecosistemicos

Mauricio Echeverry, Diego J. Lizcano, et al…

Enero, 2019