Introducción

Servicios ecosistemicos y especies amenazadas 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: Mauro

Carbono

Nutrientes

Agua

Priorización de especies

Se apilaron los modelos de distribución de Biomodelos del IAvH seleccionando dos grupos: Especies amenazadas y Mamiferos. Los modelos apilados se priorizaron espacialmente con Marxan usando portafolios de conservacion con objetivos de 10, 20, 30, 40, 50 y 75 % del habitat de las especies.

Modelos de distribución

Se descargaron 73 especies amenazadas y 116 mamiferos.

## 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

# 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


 
sp.list<-list() # crea lista vacia para almacenar especies


# all values >= -1 and <= 0.1 become 0 # fix sp chiquita problem
m <- c(-1, 0.1, 0, 0.2, 2, 1)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
# rc <- reclassify(w, rclmat)
fun1 <- function(x) { x[is.na(x)] <- 0; return(x)} # funtion to convert NA to O for Marxan


#################
# Amenazadas
#################

# Tabla atributos de las especies
tabla_sp<-read.csv("D:/BoxFiles/Box Sync/CodigoR/Servicios_Biodiversidad/data/atributos_especies.csv")


EN<-which( tabla_sp$UICN=="EN", arr.ind = TRUE)
NT<-which( tabla_sp$UICN=="NT", arr.ind = TRUE)
# LC<-which( tabla_sp$UICN=="LC", arr.ind = TRUE) # too much
VU<-which( tabla_sp$UICN=="VU", arr.ind = TRUE)
CR<-which( tabla_sp$UICN=="CR", arr.ind = TRUE)

amen <- c(EN,VU,CR)

amen_vect<-tabla_sp$scientificName[amen]
amen_list<-list()
for (i in 1:length(amen_vect)){
  sp_amen<-raster(paste("D:/BoxFiles/Box Sync/1 Analisis situacional Meta/Biodiversidad/maps/", as.character(amen_vect[i]), ".tif", sep = "" ))
  
  # rc1 <- reclassify(sp_amen, rclmat) # corrige ceros
  rc <- projectRaster(sp_amen, carbon) # reproject to put same extent
  # rm <- merge(rc, w) # extiende todos los chiquitos a w
  rc3 <- calc(rc, fun1) # remueve NAs con funcion de linea 102
  amen_list[[i]] <- rc3 # adiciona a la lista
  
}


sp_staked<- stack(amen_list) # crea el stack de la lista
names(sp_staked) <- amen_vect

# plot(sum(sp_staked), main="species richness")# grafica riqueza (suma)
sp_amen_richnes <- calc (sp_staked, fun=sum, na.rm = T)

# writeRaster(sp_amen_richnes, filename="C:/Users/diego.lizcano/Box Sync/CodigoR/Servicios_Biodiversidad/raster/stack_aenazadas.tif", format="GTiff", overwrite=TRUE)



# remueve la lista de la memoria
rm (amen_list)
rm (carbon1, carbon, nutrients, sediments, sediments1)

Unidades de conservacion y mapas de especies

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp", layer: "planning_units_wgs_84"
## with 11836 features
## It has 9 fields
## Integer64 fields read as strings:  OBJECTID_1 OBJECTID Id
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\BoxFiles\Box Sync\CodigoR\Servicios_Biodiversidad\shp\cu_servicios", layer: "cons_unit"
## with 11836 features
## It has 15 fields
cons_unit_sf <-  st_as_sf(cons_unit[10])
### costo por km2
# costo_map1<-

# Remueve rasters de memoria
# rm (cafe_rast)
# rm (papa_rast)


######################################
##### extract cost using velox raster
######################################
   
tif.vx <- velox(all.upra.geo)


extr.ls <- tif.vx$extract(sp = cons_unit_sf, fun = function(x) mean(x, na.rm = TRUE))

# Costos
# as.vector(extr.ls)



# weighted mean
# costo <- extract(all.upra, cons_unit, fun=mean)

    # # get services using the polygon hexagon map VVVERRRY SLOW !!!!!
    # # Calculates mean of service per hexagon
    # carbon_by_hexagon <- extract(carbon , cons_unit, weights=TRUE, fun=mean)
    # nutrient_by_hexagon <- extract(nutrients , cons_unit, weights=TRUE, fun=mean)
    # 

## Calculates mean of service per hexagon
# servicios_scaled <- scale(servicios) # scaled
# services_by_hexagon <- as.data.frame(extract(servicios_scaled , cons_unit, weights=TRUE, fun=mean))

# cons_unit$Carbon <- services_by_hexagon$Carbon
# cons_unit$Nutrients <- services_by_hexagon$Nutrients
# cons_unit$Sediments <- services_by_hexagon$Sediments
# cons_unit$Water <- services_by_hexagon$Water

    # 
    # 
    # # change NAS by average
    # # ind<-which(is.na(cost_by_hexagon)==TRUE)
    # # cost_by_hexagon[ind]<-960 #### put this value 1st quartil to NA
    # cons_unit$carbon<-as.vector(carbon_by_hexagon) # asign planning unit costs to hexagons
    # cons_unit$nutrient<-as.vector(nutrient_by_hexagon)
    # 
    # ##### cost uniform
    # cons_unit$cost <- 1
    # cost_by_hexagon <- cons_unit2[11]
    # 
    # # get planning unit costs from ArcGis Table
    # # costtable<-read.csv("C:/Users/Diego/Documents/CodigoR/Aldemar_anuros/shp/tablepop4.txt")
    # # cons_unit$cost<-costtable$MEAN # this is the column
    #    
    # # plot(cons_unit, col=cons_unit$cost)  ##### VERY S L O W !!!
    # # spplot(cons_unit, c("cost", "status"),edge.col = "transparent",basemap="road", grayscale=TRUE) 
    # spplot(cons_unit, c("carbon", "nutrient"), edge.col = "transparent",basemap="road", grayscale=TRUE) 
    # 
    # # delete large areas
    # cons_unit$Shape_Area <- NA
    # 
   # writeOGR(obj=cons_unit, dsn="C:/Users/diego.lizcano/Box Sync/CodigoR/Servicios_Biodiversidad/shp", layer="cu_servicios/cons_unit", driver="ESRI Shapefile") 


############## ############## ############## 
############## Mammal ############## 
############## ############## ############## 
mammal_stack <- readRDS("D:/BoxFiles/Box Sync/CodigoR/Servicios_Biodiversidad/raster/mammal_stack.rds")
mammal_stack <-  projectRaster(mammal_stack, nutrients1)
# sp_staked <- mammal_stack # change here to use mammals

mammal_stack <- calc (mammal_stack, fun=sum, na.rm = T)
names(mammal_stack) <- "Mamiferos"
mammal_stack <-  projectRaster(mammal_stack, nutrients1)
# plot cons units services
# spplot(cons_unit2, c("carbon", "nutrient"), edge.col = "transparent",basemap="road", grayscale=TRUE)

Resultados

Primer escenario, sin areas protegidas, costo uniforme

# make reserve systems 
# 
# 
cons_unit_scn1<-cons_unit


# set costs 1
cons_unit_scn1@data$cost<- 1 

# set status values
# note the 'L' after the zero is used to indicate
# that we mean the integer zero and not the decimal
# place number zero
cons_unit_scn1@data$status<-0L

# Here, we will generate a portfolio of reserve systems that represent 100%, 50% 50% of each species. 
# results<-marxan(taspu, tasinvis, targets="20%", NUMREPS=100L, BLM=0)

######### diferential targets
# Get the diferencia targets from a table IUCN
# targ<-read.csv("C:/Users/Diego/Documents/CodigoR/Aldemar_anuros/code/targets_per_sp.csv")
# target_percent<-paste(targ$target_sp,"%", sep="") 


result_scn75<-marxan(cons_unit_scn1, sp_staked, 
                     targets="75%", # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     VERBOSITY= 1L,
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2) # edge effect


sps <- amen_vect
targets75<-as.data.frame(targetsmet(result_scn75)) # Make a data frame using targets meet
sp_notargeted75<-which(apply(targets75, 2, sum)<=0) # which species do not meet any target
sp_notargete_names75<-sps[sp_notargeted75]



result_scn10<-marxan(cons_unit_scn1, sp_staked, 
                     targets= "10%",# target_percent, # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2)
targets10<-as.data.frame(targetsmet(result_scn10)) # Make a data frame using targets meet
sp_notargeted10<-which(apply(targets10, 2, sum)<=0) # which species do not meet any target
sp_notargete_names10<-sps[sp_notargeted10]



result_scn20<-marxan(cons_unit_scn1, sp_staked, 
                     targets= "20%",# target_percent, # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2)
targets20<-as.data.frame(targetsmet(result_scn20)) # Make a data frame using targets meet
sp_notargeted20<-which(apply(targets20, 2, sum)<=0) # which species do not meet any target
sp_notargete_names20<-sps[sp_notargeted20]




result_scn30<-marxan(cons_unit_scn1, sp_staked, 
                     targets= "30%",# target_percent, # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2)

targets30<-as.data.frame(targetsmet(result_scn30)) # Make a data frame using targets meet
sp_notargeted30<-which(apply(targets30, 2, sum)<=0) # which species do not meet any target
sp_notargete_names30<-sps[sp_notargeted30]


result_scn40<-marxan(cons_unit_scn1, sp_staked, 
                     targets="40%", # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2) # edge effect


targets40<-as.data.frame(targetsmet(result_scn40)) # Make a data fram using targets meet
sp_notargeted40<-which(apply(targets40, 2, sum)<=0) # which species do not meet any target
sp_notargete_names40<-sps[sp_notargeted40]



result_scn50<-marxan(cons_unit_scn1, sp_staked, 
                     targets="50%", # level of protection for each sp
                     spf=1, # species penalty factor for species
                     NUMREPS=500L, # controls the number of solutions in our portfolio
                     NCORES=3L, # number of threads for parallel processing
                     BLM=2) # edge effect


targets50<-as.data.frame(targetsmet(result_scn50)) # Make a data fram using targets meet
sp_notargeted50<-which(apply(targets50, 2, sum)<=0) # which species do not meet any target
sp_notargete_names50<-sps[sp_notargeted50]



# get levels of representation in each portfolio
results1_75.repr<-rowMeans(targetsmet(result_scn75))
results1_50.repr<-rowMeans(targetsmet(result_scn50))
results1_40.repr<-rowMeans(targetsmet(result_scn40))
results1_30.repr<-rowMeans(targetsmet(result_scn30))
results1_20.repr<-rowMeans(targetsmet(result_scn20))
results1_10.repr<-rowMeans(targetsmet(result_scn10))



# print best level of representation
# print(max(results25.repr))
print("best level of representation 75 to 10")
## [1] "best level of representation 75 to 10"
## [1] 0.5753425
## [1] 0.9178082
## [1] 0.9452055
## [1] 0.9863014
## [1] 1
## [1] 1
# histogram of proportion of species adequately
# represented in each solution
# if many of the solutions adequately represented the classes
# most of the bins would be close to 1, whereas if 
# the solutions failed to represent the classes most of the
# bins would be close to zero

# create 2 plotting areas in the one window
op<-par(mfrow=c(6,1))
  hist(rowMeans(targetsmet(result_scn10)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 10',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )
  hist(rowMeans(targetsmet(result_scn20)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 20',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )
    hist(rowMeans(targetsmet(result_scn30)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 35',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )
    hist(rowMeans(targetsmet(result_scn40)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 40',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )
        hist(rowMeans(targetsmet(result_scn50)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 50',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )
            hist(rowMeans(targetsmet(result_scn75)), freq=TRUE, xlim=c(0.6,1), las=1,
     main='Histogram of representation in portfolio 75',
     ylab='Frequency of solutions',
     xlab='Proportion of frog species adequately represented'
    )

## [1] "species 9 Pauxi pauxi in target 40%"

Especies que no cumplen con 75%: Gustavia petiolata, Juglans neotropica, Macroagelaius subalaris, Pauxi pauxi, Agamia agami, Amazona festiva, Ara militaris, Cedrela odorata, Centrolene buckleyi, Chelonoidis denticulata, Conopias cinchoneti, Conostegia superba, Elaeagia pastoensis, Grallaria alleni, Herpsilochmus axillaris, Jacaranda mimosifolia, Lagothrix lagotricha, Leptosittaca branickii, Patagioenas subvinacea, Podocnemis erythrocephala, Podocnemis unifilis, Polylepis pauta, Ramphastos ambiguus, Retrophyllum rospigliosii, Schefflera diplodactyla, Sericossypha albocristata, Tapirus terrestris, Tayassu pecari, Tinamus tao, Touit huetii, Swartzia oraria .
Especies que no cumplen con 50%: Pauxi pauxi, Cedrela odorata, Polylepis pauta, Schefflera diplodactyla, Touit huetii, Swartzia oraria.
Especies que no cumplen con 40%: Pauxi pauxi, Cedrela odorata, Polylepis pauta, Schefflera diplodactyla.
Especies que no cumplen con 30%: Schefflera diplodactyla. Especies que no cumplen con 20%: . Especies que no cumplen con 10%: .

Segundo escenario, con areas protegidas, costo uniforme

#  2. Conservacion de de 75%, 50%, 25% y mix, sin tener en cuenta el costo de oportunidad y 
# con la red de areas protegidas.
  
# geoplot richness in planning units
# with a satellite base map
# spplot(result_scn1, var='occ', basemap='satellite')


# get planning unit costs from ArcGis Table
# costtable<-read.csv("C:/Users/Diego/Documents/CodigoR/Aldemar_anuros/shp/tablepop4.txt")
# pu.costs<-costtable$MEAN # this is the column
pu.costs <- cons_unit# asign planning unit costs to hexagons as 1
pu.costs$cost <- 1
# change NAS by average
# ind<-which(is.na(pu.costs$cost)==TRUE)
# pu.costs$cost[ind]<-0 #### put 0 to NA



# get planning unit ids
pu.ids<-cons_unit@data$id
# change NAS by average
ind<-which(is.na(cons_unit_scn1@data$parque)==TRUE)
cons_unit_scn1$parque[ind]<-0 #### put 0 to NA
# change 1 by 2
is1 <- which(cons_unit_scn1@data$parque==1)
cons_unit_scn1$parque[is1]<-2 #### put 2 to 1

cons_unit_scn1$parque <- as.integer(cons_unit_scn1$parque)
# get planning unit statuses from original file
# pu.status<-cons_unit@data$status
  
# copy input parameters and data in results2, 
# change planning unit costs and statuses
# rerun MARXAN,
# load outputs into R and store them in results3
results2_75<-update(result_scn75, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque)) # cero here

results2_50<-update(result_scn50, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

results2_40<-update(result_scn40, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))


results2_30<-update(result_scn30, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

results2_20<-update(result_scn20, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))


results2_10<-update(result_scn10, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

# Lets compare 6 portfolios
# get levels of representation in each portfolio
results75.repr<-rowMeans(targetsmet(results2_75))
results50.repr<-rowMeans(targetsmet(results2_50))
results40.repr<-rowMeans(targetsmet(results2_40))
results30.repr<-rowMeans(targetsmet(results2_30))
results20.repr<-rowMeans(targetsmet(results2_20))
results10.repr<-rowMeans(targetsmet(results2_10))


targets75<-as.data.frame(targetsmet(results2_75)) # Make a data frame using targets meet
sp_notargeted75<-which(apply(targets75, 2, sum)<=0) # which species do not meet any target
sp_notargete_names75<-sps[sp_notargeted75]

targets50<-as.data.frame(targetsmet(results2_50)) # Make a data fram using targets meet
sp_notargeted50<-which(apply(targets50, 2, sum)<=0) # which species do not meet any target
sp_notargete_names50<-sps[sp_notargeted50]
  
targets40<-as.data.frame(targetsmet(results2_40)) # Make a data fram using targets meet
sp_notargeted40<-which(apply(targets40, 2, sum)<=0) # which species do not meet any target
sp_notargete_names40<-sps[sp_notargeted40] 

targets30<-as.data.frame(targetsmet(results2_30)) # Make a data fram using targets meet
sp_notargeted30<-which(apply(targets30, 2, sum)<=0) # which species do not meet any target
sp_notargete_names30<-sps[sp_notargeted30] 

targets20<-as.data.frame(targetsmet(results2_20)) # Make a data fram using targets meet
sp_notargeted20<-which(apply(targets20, 2, sum)<=0) # which species do not meet any target
sp_notargete_names20<-sps[sp_notargeted20] 

targets10<-as.data.frame(targetsmet(results2_10)) # Make a data fram using targets meet
sp_notargeted10<-which(apply(targets10, 2, sum)<=0) # which species do not meet any target
sp_notargete_names10<-sps[sp_notargeted10] 



print(paste(sp_notargete_names75, " do not meet 75% target", sep=""))
##  [1] "Gustavia petiolata do not meet 75% target"       
##  [2] "Juglans neotropica do not meet 75% target"       
##  [3] "Macroagelaius subalaris do not meet 75% target"  
##  [4] "Pauxi pauxi do not meet 75% target"              
##  [5] "Agamia agami do not meet 75% target"             
##  [6] "Amazona festiva do not meet 75% target"          
##  [7] "Ara militaris do not meet 75% target"            
##  [8] "Cedrela odorata do not meet 75% target"          
##  [9] "Centrolene buckleyi do not meet 75% target"      
## [10] "Chelonoidis denticulata do not meet 75% target"  
## [11] "Conopias cinchoneti do not meet 75% target"      
## [12] "Conostegia superba do not meet 75% target"       
## [13] "Elaeagia pastoensis do not meet 75% target"      
## [14] "Grallaria alleni do not meet 75% target"         
## [15] "Herpsilochmus axillaris do not meet 75% target"  
## [16] "Jacaranda mimosifolia do not meet 75% target"    
## [17] "Lagothrix lagotricha do not meet 75% target"     
## [18] "Leptosittaca branickii do not meet 75% target"   
## [19] "Patagioenas subvinacea do not meet 75% target"   
## [20] "Podocnemis erythrocephala do not meet 75% target"
## [21] "Podocnemis unifilis do not meet 75% target"      
## [22] "Polylepis pauta do not meet 75% target"          
## [23] "Ramphastos ambiguus do not meet 75% target"      
## [24] "Retrophyllum rospigliosii do not meet 75% target"
## [25] "Schefflera diplodactyla do not meet 75% target"  
## [26] "Sericossypha albocristata do not meet 75% target"
## [27] "Tapirus terrestris do not meet 75% target"       
## [28] "Tayassu pecari do not meet 75% target"           
## [29] "Tinamus tao do not meet 75% target"              
## [30] "Touit huetii do not meet 75% target"             
## [31] "Swartzia oraria do not meet 75% target"
## [1] "Pauxi pauxi do not meet 50% target"            
## [2] "Cedrela odorata do not meet 50% target"        
## [3] "Polylepis pauta do not meet 50% target"        
## [4] "Schefflera diplodactyla do not meet 50% target"
## [5] "Touit huetii do not meet 50% target"           
## [6] "Swartzia oraria do not meet 50% target"
## [1] "Pauxi pauxi do not meet 40% target"            
## [2] "Cedrela odorata do not meet 40% target"        
## [3] "Polylepis pauta do not meet 40% target"        
## [4] "Schefflera diplodactyla do not meet 40% target"
## [1] "Schefflera diplodactyla do not meet 30% target"
## [1] " do not meet 20% target"
## [1] " do not meet 10% target"
# create 4 plotting areas in the one window
op2<-par(mfrow=c(6,1))

  # histogram of first portfolio
  hist(results75.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 75%"
  )
  
  # print best level of representation
  # print(max(results100.repr))
  
  # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results50.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 50%"
  )
  
  # print best level of representation
  # print(max(results50.repr))
  
    # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results40.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 40%"
  )
  
    hist(results30.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 30%"
  )
  
  hist(results20.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 20%"
  )
  
    hist(results10.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 10%"
  )

## [1] "best level of representation"
## [1] 0.5753425
## [1] 0.9178082
## [1] 0.9452055
## [1] 0.9863014
## [1] 1
## [1] 1
## [1] "best solution"

## [1] "species 9 Pauxi pauxi in target 40%"

## [1] "results2_10 not shown, for some reason produces error"

# results2_75_sf <- st_as_sf(results2_75)
  # plot(results2_10, colramp='YlGnBu')#, basemap="road", grayscale=TRUE, force_reset = FALSE)
# The solutions in this portfolio are fairly fragmented. 
# All the solutions so far were made under the assumption that all planning units have equal acquisition costs 
# and that Col does not have any protected areas. Not true!
  
  
##### ##### ##### ##### 
##### save files figs
##### ##### ##### ##### 
# 1. Open jpeg file
#  jpeg("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.jpg", width = 2500, height = 1750)
# 
#  png("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.png", width = 2500, height = 1760, bg="white")
# 
#  pdf("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.pdf", width = 70, height = 50)
# # 2. Create the plot
#  plot(results3_30)#, 0)
# # 3. Close the file
#  dev.off()


# # #### write result
#  write.MarxanSolved(results3_10, "E:/Orinoquia/E3/10")
#  write.MarxanSolved(results3_20, "E:/Orinoquia/E3/20")
#  write.MarxanSolved(results3_30, "E:/Orinoquia/E3/30")
#  write.MarxanSolved(results3_40, "E:/Orinoquia/E3/40")
#  write.MarxanSolved(results3_50, "E:/Orinoquia/E3/50")
#  write.MarxanSolved(results3_75, "E:/Orinoquia/E3/75")
# #   

# ########################
# ## Convert to shp ####
# ########################
# 
# selection2_75 <- (t(results2_75@results@selections)) # get selection
# selection2_75 <- as.data.frame(1*selection2_75) # convert to number
# names(selection2_75) <- 1:dim(selection2_75)[[2]] # put names
# row.names (selection2_75) <- 1:dim(selection2_75)[[1]] # put names
# selection2_75$id <- 1:dim(selection2_75)[[1]]
# selection2_75$frec_sel <- apply(selection2_75,1,mean) # sum
# cons_unit_sf2 <-  merge(cons_unit_sf, selection2_75,
#                         by.x = "id", by.y = "id") # paste to sf
# 
# st_write(cons_unit_sf2, "D:/BoxFiles/Box Sync/CodigoR/Servicios_Biodiversidad/shp/E2/75_selection.shp" , driver = "ESRI Shapefile")

Tercer escenario, con areas protegidas y costo UPRA

#  2. Conservacion de de 75%, 50%, 25% y mix, sin tener en cuenta el costo de oportunidad y 
# con la red de areas protegidas.
  
# geoplot richness in planning units
# with a satellite base map
# spplot(result_scn1, var='occ', basemap='satellite')
# make reserve systems 
# 
# 
cons_unit_scn1<-cons_unit


# set costs 1

# get planning unit costs from ArcGis Table
# costtable<-read.csv("C:/Users/Diego/Documents/CodigoR/Aldemar_anuros/shp/tablepop4.txt")
# pu.costs<-costtable$MEAN # this is the column
pu.costs <- cons_unit# asign planning units

# Using cost extracted from UPRA
pu.costs$cost <- as.vector(extr.ls)
# change NAS by average
ind<-which(is.na(pu.costs$cost)==TRUE)
pu.costs$cost[ind]<-0 #mean(as.vector(extr.ls), na.rm = TRUE) #### put mean cost to NA



# get planning unit ids
pu.ids<-cons_unit@data$id
# change NAS by average
ind<-which(is.na(cons_unit_scn1@data$parque)==TRUE)
cons_unit_scn1$parque[ind]<-0 #### put 0 to NA
# change 1 by 2
is1 <- which(cons_unit_scn1@data$parque==1)
cons_unit_scn1$parque[is1]<-2 #### put 2 to 1

cons_unit_scn1$parque <- as.integer(cons_unit_scn1$parque)
# get planning unit statuses from original file
# pu.status<-cons_unit@data$status
  


# copy input parameters and data in results2, 
# change planning unit costs and statuses
# rerun MARXAN,
# load outputs into R and store them in results3
results3_75<-update(result_scn75, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque)) # cero here

results3_50<-update(result_scn50, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

results3_40<-update(result_scn40, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))


results3_30<-update(result_scn30, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

results3_20<-update(result_scn20, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))


results3_10<-update(result_scn10, ~pu(pu.ids, 
                                        cost=pu.costs$cost, 
                                        status=cons_unit_scn1$parque))

# Lets compare 6 portfolios
# get levels of representation in each portfolio
results75.repr<-rowMeans(targetsmet(results3_75))
results50.repr<-rowMeans(targetsmet(results3_50))
results40.repr<-rowMeans(targetsmet(results3_40))
results30.repr<-rowMeans(targetsmet(results3_30))
results20.repr<-rowMeans(targetsmet(results3_20))
results10.repr<-rowMeans(targetsmet(results3_10))


targets75<-as.data.frame(targetsmet(results3_75)) # Make a data frame using targets meet
sp_notargeted75<-which(apply(targets75, 2, sum)<=0) # which species do not meet any target
sp_notargete_names75<-sps[sp_notargeted75]

targets50<-as.data.frame(targetsmet(results3_50)) # Make a data fram using targets meet
sp_notargeted50<-which(apply(targets50, 2, sum)<=0) # which species do not meet any target
sp_notargete_names50<-sps[sp_notargeted50]
  
targets40<-as.data.frame(targetsmet(results3_40)) # Make a data fram using targets meet
sp_notargeted40<-which(apply(targets40, 2, sum)<=0) # which species do not meet any target
sp_notargete_names40<-sps[sp_notargeted40] 

targets30<-as.data.frame(targetsmet(results3_30)) # Make a data fram using targets meet
sp_notargeted30<-which(apply(targets30, 2, sum)<=0) # which species do not meet any target
sp_notargete_names30<-sps[sp_notargeted30] 

targets20<-as.data.frame(targetsmet(results3_20)) # Make a data fram using targets meet
sp_notargeted20<-which(apply(targets20, 2, sum)<=0) # which species do not meet any target
sp_notargete_names20<-sps[sp_notargeted20] 

targets10<-as.data.frame(targetsmet(results3_10)) # Make a data fram using targets meet
sp_notargeted10<-which(apply(targets10, 2, sum)<=0) # which species do not meet any target
sp_notargete_names10<-sps[sp_notargeted10] 



print(paste(sp_notargete_names75, " do not meet 75% target", sep=""))
##  [1] "Gustavia petiolata do not meet 75% target"       
##  [2] "Juglans neotropica do not meet 75% target"       
##  [3] "Macroagelaius subalaris do not meet 75% target"  
##  [4] "Pauxi pauxi do not meet 75% target"              
##  [5] "Agamia agami do not meet 75% target"             
##  [6] "Amazona festiva do not meet 75% target"          
##  [7] "Ara militaris do not meet 75% target"            
##  [8] "Cedrela odorata do not meet 75% target"          
##  [9] "Centrolene buckleyi do not meet 75% target"      
## [10] "Chelonoidis denticulata do not meet 75% target"  
## [11] "Conopias cinchoneti do not meet 75% target"      
## [12] "Conostegia superba do not meet 75% target"       
## [13] "Elaeagia pastoensis do not meet 75% target"      
## [14] "Grallaria alleni do not meet 75% target"         
## [15] "Herpsilochmus axillaris do not meet 75% target"  
## [16] "Jacaranda mimosifolia do not meet 75% target"    
## [17] "Lagothrix lagotricha do not meet 75% target"     
## [18] "Leptosittaca branickii do not meet 75% target"   
## [19] "Patagioenas subvinacea do not meet 75% target"   
## [20] "Podocnemis erythrocephala do not meet 75% target"
## [21] "Podocnemis unifilis do not meet 75% target"      
## [22] "Polylepis pauta do not meet 75% target"          
## [23] "Ramphastos ambiguus do not meet 75% target"      
## [24] "Retrophyllum rospigliosii do not meet 75% target"
## [25] "Schefflera diplodactyla do not meet 75% target"  
## [26] "Sericossypha albocristata do not meet 75% target"
## [27] "Tapirus terrestris do not meet 75% target"       
## [28] "Tayassu pecari do not meet 75% target"           
## [29] "Tinamus tao do not meet 75% target"              
## [30] "Touit huetii do not meet 75% target"             
## [31] "Swartzia oraria do not meet 75% target"
## [1] "Pauxi pauxi do not meet 50% target"            
## [2] "Agamia agami do not meet 50% target"           
## [3] "Ara militaris do not meet 50% target"          
## [4] "Cedrela odorata do not meet 50% target"        
## [5] "Polylepis pauta do not meet 50% target"        
## [6] "Schefflera diplodactyla do not meet 50% target"
## [7] "Touit huetii do not meet 50% target"           
## [8] "Swartzia oraria do not meet 50% target"
## [1] "Pauxi pauxi do not meet 40% target"            
## [2] "Agamia agami do not meet 40% target"           
## [3] "Cedrela odorata do not meet 40% target"        
## [4] "Polylepis pauta do not meet 40% target"        
## [5] "Schefflera diplodactyla do not meet 40% target"
## [1] "Cedrela odorata do not meet 30% target"        
## [2] "Schefflera diplodactyla do not meet 30% target"
## [1] "Schefflera diplodactyla do not meet 20% target"
## [1] " do not meet 10% target"
# create 4 plotting areas in the one window
op2<-par(mfrow=c(6,1))

  # histogram of first portfolio
  hist(results75.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 75%"
  )
  
  # print best level of representation
  # print(max(results100.repr))
  
  # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results50.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 50%"
  )
  
  # print best level of representation
  # print(max(results50.repr))
  
    # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results40.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 40%"
  )
  
    hist(results30.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 30%"
  )
  
  hist(results20.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 20%"
  )
  
    hist(results10.repr, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 10%"
  )

## [1] "best level of representation"
## [1] 0.5753425
## [1] 0.890411
## [1] 0.9315068
## [1] 0.9726027
## [1] 0.9863014
## [1] 1
## [1] "best solution"

## [1] "species 9 Pauxi pauxi in target 40%"

##### ##### ##### ##### 
##### save files figs
##### ##### ##### ##### 
# 1. Open jpeg file
#  jpeg("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.jpg", width = 2500, height = 1750)
# 
#  png("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.png", width = 2500, height = 1760, bg="white")
# 
#  pdf("E:/Orinoquia/E3/fig/result_scn_sel_frec_30.pdf", width = 70, height = 50)
# # 2. Create the plot
#  plot(results3_30)#, 0)
# # 3. Close the file
#  dev.off()


# # #### write result
#  write.MarxanSolved(results3_10, "E:/Orinoquia/E3/10")
#  write.MarxanSolved(results3_20, "E:/Orinoquia/E3/20")
#  write.MarxanSolved(results3_30, "E:/Orinoquia/E3/30")
#  write.MarxanSolved(results3_40, "E:/Orinoquia/E3/40")
#  write.MarxanSolved(results3_50, "E:/Orinoquia/E3/50")
#  write.MarxanSolved(results3_75, "E:/Orinoquia/E3/75")
# #   

# ########################
# ## Convert to shp ####
# ########################
# 
# selection2_10 <- (t(results3_10@results@selections)) # get selection
# selection3_10 <- as.data.frame(1*selection3_10) # convert to number
# names(selection3_10) <- 1:dim(selection3_10)[[2]] # put names
# row.names (selection3_10) <- 1:dim(selection3_10)[[1]] # put names
# selection3_10$id <- 1:dim(selection3_10)[[1]]
# selection3_10$frec_sel <- apply(selection3_10,1,mean) # sum
# cons_unit_sf2 <-  merge(cons_unit_sf, selection3_10,
#                         by.x = "id", by.y = "id") # paste to sf
# 
# st_write(cons_unit_sf2, "D:/BoxFiles/Box Sync/CodigoR/Servicios_Biodiversidad/shp/E3/10_selection.shp" , driver = "ESRI Shapefile")
#

Cuarto escenario, con areas protegidas, costo UPRA y target (%) variable, dependiendo de la especie. Aun Pendiente

#  3. Conservacion de de 75%, 50%, 25% y sin teniendo en cuenta el costo de oportunidad y 
# teniendo en cuenta red actual de areas protegidas.
  
# geoplot richness in planning units
# with a satellite base map
# spplot(result_scn1, var='occ', basemap='satellite')


# get planning unit costs from ArcGis Table
# costtable<-read.csv("C:/Users/Diego/Documents/CodigoR/Aldemar_anuros/shp/tablepop4.txt")
# pu.costs<-as.numeric(costtable$MEAN) # this is the column

pu.costs<-cost_by_hexagon


# get planning unit ids
pu.ids<-cons_unit@data$id 
PAs<-as.integer(cons_unit2$parque) # get protected areas 
# get planning unit statuses from original file
# pu.status<-cons_unit@data$status
  
# copy input parameters and data in results2, 
# change planning unit costs and statuses
# rerun MARXAN,
# load outputs into R and store them in results3
results3_75<-update(results2_75, ~pu(id=pu.ids, 
                                        cost=pu.costs, 
                                        status=PAs)) 

results3_50<-update(results2_50, ~pu(id=pu.ids, 
                                        cost=pu.costs, 
                                        status=PAs))

results3_25<-update(results2_25, ~pu(id=pu.ids, 
                                        cost=pu.costs, 
                                        status=PAs))

results3_mix<-update(results2_mix, ~pu(id=pu.ids, 
                                        cost=pu.costs, 
                                        status=PAs))



# Lets compare 3 portfolios
# get levels of representation in each portfolio
results75.repr3<-rowMeans(targetsmet(results3_75))
results50.repr3<-rowMeans(targetsmet(results3_50))
results25.repr3<-rowMeans(targetsmet(results3_25))
resultsmix.repr3<-rowMeans(targetsmet(results3_mix))



targets75<-as.data.frame(targetsmet(results3_75)) # Make a data frame using targets meet
sp_notargeted75<-which(apply(targets75, 2, sum)<=0) # which species do not meet any target
sp_notargete_names75<-sps[sp_notargeted75]

targets50<-as.data.frame(targetsmet(results3_50)) # Make a data fram using targets meet
sp_notargeted50<-which(apply(targets50, 2, sum)<=0) # which species do not meet any target
sp_notargete_names50<-sps[sp_notargeted50]
  
targets25<-as.data.frame(targetsmet(results3_25)) # Make a data fram using targets meet
sp_notargeted25<-which(apply(targets25, 2, sum)<=0) # which species do not meet any target
sp_notargete_names25<-sps[sp_notargeted25] 

targetsdif<-as.data.frame(targetsmet(results3_mix)) # Make a data frame using targets meet
sp_notargeteddif<-which(apply(targetsdif, 2, sum)<=0) # which species do not meet any target
sp_notargete_namesdif<-sps[sp_notargeteddif]

print(paste(sp_notargete_names75, " do not meet 75% target", sep=""))
print(paste(sp_notargete_names50, " do not meet 50% target", sep=""))
print(paste(sp_notargete_names25, " do not meet 25% target", sep=""))
print(paste(sp_notargete_namesdif, " do not meet mix% target", sep=""))



# create 3 plotting areas in the one window
op3<-par(mfrow=c(4,1))

  # histogram of first portfolio
  hist(results75.repr3, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 75%"
  )
  
  # print best level of representation
  #print(max(results100.repr))
  
  # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results50.repr3, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 50%"
  )
  
  # print best level of representation
  # print(max(results50.repr))
  
    # histogram of second portfolio
  # if you see a giant single rectangle this means
  # all the solutions have the same level of representation
  hist(results25.repr3, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with 25%"
  )
    hist(resultsmix.repr3, freq=TRUE, xlim=c(0,1), las=1,
       ylab='Frequency of solutions',
       xlab='Proportion of species adequately represented',
       main="Level of representation with mix%"
  )

 par(op2) 
  
    
  # print best level of representation
  # print(max(results25.repr))
  

  # make a geoplot of the best solution
  plot(results3_75, 0, basemap="road", grayscale=TRUE, force_reset = FALSE)
  plot(results3_50, 0, basemap="road", grayscale=TRUE, force_reset = FALSE)
  plot(results3_25, 0, basemap="road", grayscale=TRUE, force_reset = FALSE)
  plot(results3_mix, 0, basemap="road", grayscale=TRUE, force_reset = FALSE)
  
  # planing uniit selection frecuency
  plot(results3_75, colramp='YlGnBu', basemap="road", grayscale=TRUE, force_reset = FALSE)
  plot(results3_50, colramp='YlGnBu', basemap="road", grayscale=TRUE, force_reset = FALSE)
  plot(results3_25, colramp='YlGnBu', basemap="road", grayscale=TRUE, force_reset = FALSE)  
  plot(results3_mix, colramp='YlGnBu', basemap="road", grayscale=TRUE, force_reset = FALSE) 
  
# The solutions in this portfolio are fairly fragmented. 
# All the solutions so far were made under the assumption that all planning units have equal acquisition costs 
# and that Col does not have any protected areas. Not true!

  
targets75<-as.data.frame(targetsmet(results3_75)) # Make a data frame using targets meet
sp_notargeted75<-which(apply(targets75, 2, sum)<=0) # which species do not meet any target
sp_notargete_names75<-sps[sp_notargeted75]

targets50<-as.data.frame(targetsmet(results3_50)) # Make a data fram using targets meet
sp_notargeted50<-which(apply(targets50, 2, sum)<=0) # which species do not meet any target
sp_notargete_names50<-sps[sp_notargeted50]
  
targets25<-as.data.frame(targetsmet(results3_25)) # Make a data fram using targets meet
sp_notargeted25<-which(apply(targets25, 2, sum)<=0) # which species do not meet any target
sp_notargete_names25<-sps[sp_notargeted25] 

targetsdif<-as.data.frame(targetsmet(results3_mix)) # Make a data frame using targets meet
sp_notargeteddif<-which(apply(targetsdif, 2, sum)<=0) # which species do not meet any target
sp_notargete_namesdif<-sps[sp_notargeteddif]

print(paste(sp_notargete_names75, " do not meet 75% target", sep=""))
print(paste(sp_notargete_names50, " do not meet 50% target", sep=""))
print(paste(sp_notargete_names25, " do not meet 25% target", sep=""))
print(paste(sp_notargete_namesdif, " do not meet mix% target", sep=""))

Escenario cinco variando el BLM, para E3 target 20%

# All the solutions in our current portfolio, results3, seem to be fairly fragmented. 
# If implemented as protected areas, these solutions might be associated with poor connectivity 
# and high maintenance costs. To reduce fragmentation, we can increase the boundary length
# modifier (BLM). However, in order to maintain adequate levels of representation for the 
# species, MARXAN will select more expensive planning units. 
# How can we pick an appropriate BLM while still making sure the acquisition costs are adequate cost? 
# Let's generate six more portfolios, each using a different BLM, and plot the trade-off 
# between acquisition cost and fragmentation using the best solutions in each portfolio.

## generate list of portfolios with different BLMS
# make vector BLM parameters to use
blm.pars=c(1, 10, 100, 250, 500, 750, 1000)


# create list with different portfolio for each BLM for results2_10
results4_mix<-list()
for (i in seq_along(blm.pars)) {
  results4_mix[[i]]<-update(results3_20, ~opt(BLM=blm.pars[i], NUMREPS=10L))
}

## extract data from portfolios
# create empty vectors to store values
cost<-c()
con<-c()
blm<-c()

# extract values for best solutions
for (i in seq_along(blm.pars)) {
  cost<-append(cost, summary(results4_mix[[i]])[["Cost"]])
  con<-append(con, summary(results4_mix[[i]])[["Connectivity"]])
  blm<-append(blm, rep(blm.pars[i], nrow(summary(results4_mix[[i]]))))
}

## plot trade-off between shortfall and connectivity
# get colours for legend
legend.cols<-c("#ffffe2", "#FFFFB2", "#FED976", "#FEB24C", "#FD8D3C", "#F03B20", "#BD0026")
pt.cols<-legend.cols[match(blm, blm.pars)]

# reset plotting window
# par(mfrow=c(1,1))

# plot trade-off data
# higher shortfall values means worse representation
# higher connectivity values mean more fragmentation
plot(cost~con, bg=pt.cols, col='black', ylab='Cost', xlab='Connectivity', pch=21,
     main='Trade-off between cost and connectivity')
abline(lm(cost~con))

# add legend
legend("topright", legend=blm.pars, col='black', pt.bg=legend.cols, pch=21, title='BLM')