Análisis de las causas de transformación de los ecosistemas naturales en el departamento del Meta 2000-2014

INTRODUCCIÓN

1. Transformación antrópica de ecosistemas naturales

Los patrones historicos de la distribución espacial y la extensión de los ecosistemas en Colombia son el producto de la interacción entre los ecosistemas naturales y el proceso de poblamiento historico del país (Etter et al. 2008). En las decadas recientes los cambios en las coberturas y usos de la tierra se han acelerado producto del incremento poblacional y los niveles de consumo, que se han convertido en los mayores factores de cambio global (Vitousek et al. 1997 Foley et al. (2005)). Usualmente, los cambios en coberturas se plantea como un proceso de incremento lineal en la transformación de condiciones naturales a transformadas (Lambin et al. 2003, Sánchez-Cuervo et al. (2012)). Sin embargo, la transformación antrópica de los ecosistemas naturales pueden tener periodos donde las dinámicas cambian de acuerdo a los cambios poblacionales (Armenteras et al. 2011).

La colonización española desde los 1500 representa el proceso que a más largo plazo dio inicio a la transformación de ecosistemas en Colombia, donde inicialmente la expansión de tierras para la ganadería y recientemente la ampliación de la frontera agrícola son actividades que han tenido mayor impacto en los ecosistemas naturales a largo plazo (Etter et al. 2008). Despues de los 1700, la huella espacial de la agricultura se aumento continuamente con el incremento de la población humana y de manera acelerada en los 1900, lo que resulto en la pérdida de un tercio de la extensión total de los ecosistemas forestales natrurales de Colombia, donde la actividad de ganadería jugo el papel principal en el proceso de transformación (Armenteras et al. 2013a). Sin embargo, ha habido una sucesión en los factores historicos como es el caso entre 1850 y 1920 donde la apertura economica a mercados internacionales estimulo la expansión de la actividad agrícola y la extracción de recursos naturales lo que introdujo factores adicionales a la ganadería (Armenteras et al. 2013b).

La pérdida de los bosques tienen efectos directos obre la pérdida de biodiversidad (Brooks et al. 2002 Jha and Bawa (2006)). En Colombia las causas próximas varían por región y han sido descritas por (Armenteras et al. 2013c) en las subregiones: piedemonte, norte, nororiente y sur. El piedemonte es la zona de más desarrollo económico de la región y concentra la mayor parte de la deforestación debida principalmente a la colonización y a actividades de ganadería, minería y cultivos ilícitos. La subregión norte carece de alternativas económicas viables y allí las zonas de reserva forestal han sido ocupadas con la expectativa de sustracción y titulación, fortalecidas por la perspectiva de yacimientos minero-energéticos. Según el Sistema de Monitoreo de Bosques y Carbono de Colombia (SMByC) el departamento del Meta para el año 2016 estima 3.077.470 hectareas de bosque y se han pérdido 853.103 hectareas entre 1990 y 2016.

Pese a que los ecosistemas de sabanas naturales representan un cuarto de la extensión de los ecosistemas naturales de Colombia estos no han tenido la misma atención en cuanto a su conservación comparado con ecosistemas mejor estudiados como los ecosistemas de bosques andinos y amazonicos (Furley 1999). Los ecosistemas de sabanas de los llanos orientales representan el segundo sistema de sabanas más extenso despues del cerrado (Romero et al. 2009) y son de importancia internacional al contener más del 55% de los humedales del país, 40% de las aguas subterráneas del país (IDEAM 2010), 46% de los peces colombianos (Maldonado-Ocampo et al. 2013) y 40% de aves (Acevedo-Charry et al. 2014). Además, forman importantes corredores para muchas especies de mamíferos y reptiles, como el jaguar, que migra entre las regiones de Amazonas, Guayana y Andes (Boron et al. 2016). Para la los llanos orientales de Colombia se han pérdido 211.600 hectareas de ecosistemas de sabana entre los años 1987 y 2007 donde una fracción importante de los cambios de coberturas/usos estaba asociada a la expansión de cultivos de palma (Romero-Ruiz et al. 2012).

Para poder diseñar e implementar estrategias y programas efectivos y eficientes para reducir la degradación y pérdida de ecosistemas naturales es necesario entender sus causas, ya que los contextos socioeconómicos, institucionales, legales y ambientales varian regional, nacional y localmente. En el caso de ecosistemas de bosque se han cuestionado las estrategias para reducir la deforestación y degradación de bosques debido la falta de salvaguardas, de enfoque en cuanto a los beneficios de los programas y análisis de causas de deforestación.

El presente trabajo tiene como objetivo realizar un análisis para explorar las causas recientes de pérdida de ecosistemas de sabana y bosques en el departamento del Meta (2000-2014), analizando los patrones espaciales y temporales de transformación de ecosistemas de bosque y sabanas respecto a la variación de factores sociales, económicos, políticos y culturales. Así como desarrollar la base de datos socioecoómica para caracterizar los municipios del Meta.

2. Enfoque del análisis de causas tranformación de ecosistemas naturales en Meta

Para este analisis se siguio el marco conceptual propuesto por (Geist and Lambin 2001), agrupando las variables de causas próximas y subyacentes por factores poblacionales, económicos, infraestructura, tecnológicos y políticos Figura 1. (Geist and Lambin 2001) hacen una revisión de la literatura sobre causas y agentes de la deforestación en el trópico, donde las causas próximas son las actividades humanas que afectan directamente el bosque, como por ejemplo: la agricultura, la tala, la minería, y la expansión de infraestructura, mientras que las causas subyacentes se entienden como procesos fundamentales que afectan las causas próximas y que pueden ser de origen social, político, económico, y cultural. Este trabajo explora las causas próximas y subyacentes de la transformación de ecosistemas de bosques y sabanas en el departamento del Meta.

La mayoria de los estudios de las causas de la transformación de ecosistemas naturales se han enfocado en el estudio de la deforestación y degradación de bosques. Debido a que los ecosistemas de sabanas son resperesentativos dentro de los ecosistemas naturales en el departamento del Meta se realizó un análisis que se enfoco a la transformación de los ecosistemas de sabanas y ecosistemas de bosque. Para cuantificar la pérdida de bosques se utilizaron los datos generados por Hansen et al. (2013) que estiman las áreas de deforestación anual a partir de imágenes del sensor LandSat entre 2000-2014. Para estimar la transformación de sabanas se utilizaron los datos de sitios activos de fuego a generados a partir del sensor MODIS (Giglio et al. 2003) que corresponden a los sitios detectados como fuego en un área de 1km2 donde se toman medidas cuatro veces al día. Se seleccionarón los fuegos como proxy de la transformación de ecosistemas de sabana ya que el fenomeno de fuegos esta usualmente asociado a actividades antrópicas (80%) y raramente es un proceso natural (20%), como es el caso de quemas para producir pastos para la gandería.

Para realizar el análisis de causas de transformación de ecosistemas de bosques y sabanas en el departamento del Meta se utilizaron los municipios como unidades de análisis para las que se estimaron las áreas anuales de deforestación y el número de eventos anuales de fuegos. Se definiio el periodo de análisis de 2000-2014 para explorar las causas recientes de transformación y complementar análisis previos para la orinoquía tales como (Romero-Ruiz et al. 2012). Adicionalmente, como parte de la fase de revisión y selección de las variables explicativas asociadas a los diferentes factores propuesto por (Geist and Lambin 2001), se realizó la caracterización de los municipios del departamento del Meta respecto a las variables exploradas.

Diagrama de flujo producto de la sistematización y generalización de los patrones causantes de la deforestación en el trópico para 152 estudios tomado de Geist and Lambin (2001)

Diagrama de flujo producto de la sistematización y generalización de los patrones causantes de la deforestación en el trópico para 152 estudios tomado de Geist and Lambin (2001)

METODOLOGÍA

La figura 2 describe el flujo de análisis para explorar la correlación de variables socio-económicas con la transformación de ecosistemas de bosques y sabanas en el Meta; y desarrollar la base de datos para la caracterización socio-económica del departamento. A continuación se describen las variables utilizadas y los procesos de análisis aplicados.

library (sp)  # classes for spatial data
library (raster)  # grids, rasters
library (rasterVis)
library (rasterVis)  # raster visualisation
library (maptools) # map simple operations
library (rgeos) # coord transform
library (rgdal) # coord transform
library (ggmap) # visualization
library (tidyverse) # date managemgemt
library (mapview) # Visualization widget in esri sat image and open street map
library (readxl)# read excel data
library (sf) # shp and vector map handling 
library (gfcanalysis) # Hansen et al. 2013 Global Forest Change dataset
library (rgdal)# used to read/write shapefiles and rasters
library (dplyr)
library (printr)
library (corrplot) # correlations
library (knitr)
library (kableExtra)
##################################
### Read Data set from excel
##################################

# Specify sheet with a number or name
info_meta <- read_excel ("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/NASCA000228-2017_ProductoB_municipal_dic2017.xlsx", sheet = "BASE_DATOS")

meta_datos <- read_excel ("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/NASCA000228-2017_ProductoB_municipal_dic2017.xlsx", sheet = "VARIABLES")

x <- unique(meta_datos[4])
# knitr::kable(x, digits = 2, caption = "Tipos de variables y grupos")

# kable(info_meta, "html") %>%
#  kable_styling() %>%
#  scroll_box(width = "650px", height = "350px")

# kable(head(info_meta))#, "html") %>%
  #kable_styling(bootstrap_options = c("striped","hover","condensed","responsive", full_width = F)) %>%
#scroll_box(width = "650px", height = "350px")

Deforestación

Col_full<-getData('GADM', country='COL', level=2)

# Opción 2 #utm zona 17 sur con datum WGS84 
utm18 <- "+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"


# Indicate where we want to save GFC tiles downloaded from Google. For any 
# given AOI, the script will first check to see if these tiles are available 
# locally (in the below folder) before downloading them from the server - so I 
# recommend storing ALL of your GFC tiles in the same folder. For this example 
# we will save files in the current working directory folder.
if (require(snow)) beginCluster()

data_folder <- 'C:/Users/diego.lizcano/Box Sync/CodigoR/Toshiba/Amazon/HansenData'

###############################################################################
# Download data from Google server for a given AOI
###############################################################################
col_full_sf <- st_as_sf(Col_full) # convert to sf object

meta_municip <- filter(col_full_sf, NAME_1 == "Meta")

ggplot(meta_municip) +
        geom_sf(aes(fill = NAME_2))

# col_full_sf %>% filter(NAME_1 == "Meta") %>% plot()

deforelist<-list()
for (i in 1:length(unique(meta_municip$NAME_2))){ # loop statrt  
  meta_adm2_name<-unique(meta_municip$NAME_2)[i] # get name from i
  meta_adm2_poly<-meta_municip[i,] # get polygon from i
  meta_adm2_poly_utm<-st_transform(meta_adm2_poly, utm18)# transformación UTM
  
  
  # Load a demo AOI from the P drive - notice that first we specify the folder 
  # the shapefile is in, and then the name of the shapefile without the '.shp'
  # aoi <- readOGR(system.file('extdata', package='gfcanalysis'), 'ZOI_NAK_2012')
  aoi<- as(meta_adm2_poly, "Spatial")
  
  # Calculate the google server URLs for the tiles needed to cover the AOI
  tiles <- calc_gfc_tiles(aoi)
  
  # Check to see if these tiles are already present locally, and download them if 
  # they are not.
  download_tiles(tiles, data_folder)
  
  # Extract the GFC data for this AOI from the downloaded GFC tiles, mosaicing 
  # multiple tiles as necessary (if needed to cover the AOI), and saving  the 
  # output data to a GeoTIFF (can also save in ENVI format, Erdas format, etc.).
  gfc_data <- extract_gfc(aoi, data_folder, 
                          filename=paste("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/meta/",meta_adm2_name, "_extract.tif",sep = ""),
                          overwrite=TRUE)
  
  ###############################################################################
  # Performing thresholding and calculate basic statistics
  ###############################################################################
  
  # Calculate and save a thresholded version of the GFC product
  gfc_thresholded <- threshold_gfc(gfc_data, forest_threshold=90, # 90 default
                                   filename=paste("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/meta/",meta_adm2_name, "_thresholded.tif", sep = ""),
                                   overwrite=TRUE)
  
  gfc_thresholded_utm<-projectRaster(gfc_thresholded, crs=utm18) 
  
  
  
  # Calculate annual statistics on forest loss/gain
  gfc_stats <- gfc_stats(aoi, gfc_thresholded_utm)
  deforelist[[i]]<-list(meta_adm2_name, gfc_stats)
  
  # Save statistics to CSV files for use in Excel, etc.
  write.csv(gfc_stats$loss_table, file=paste("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/meta/table/",meta_adm2_name, "_losstable.csv",sep = ""), row.names=FALSE)
  write.csv(gfc_stats$gain_table, file=paste("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/meta/table/",meta_adm2_name, "_gaintable.csv",sep = ""), row.names=FALSE)
  
  
} # end loop


# Stop the parallel processing cluster
if (require(snow)) endCluster()

# embeed in a simpler list
muni_list<-list()
for (i in 1: length(deforelist)){
  depto_table<-deforelist[[i]][[2]][[1]]
  depto_table$depto<-deforelist[[i]][[1]]
  muni_list[[i]]<-depto_table
}

library(plyr)
depto_table<-ldply(muni_list)
write.csv(depto_table, "C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/municipi_table.csv")

Los datos de deforestación anual entre el periodo 2010-2014 para el departamento del Meta corresponden a los datos del Global Forest Change Dataset Hansen et al. (2013), se recortaron para cada municipio en cada año para conformar una serie de datos de tipo panel Jackson (2011). Las datos anuales de deforestación se obtuvieron a través del paquete gfcanalysis https://github.com/azvoleff/gfcanalysis con una resolución espacial de 30*30 metros. Se evaluo el uso de los datos de cuantificación de boques producto de la fuente de datos del Sistema de Monitoreo de Bosques y Carbono (SMByC) pero se opto utilizar los datos de Hansen et al. (2013) ya que son datos anuales durante el periodo 2000-2014 y los datos del SMByC corresponden tienen vacios de información anual para los años: 2001, 2002, 2003, 2004, 2006, 2007, 2008 y 2009.

A. Sitios detectados de quemas anuales (Giglio et al. 2003) para el departamento del Meta y B. Extensión de áreas deforestadas para el departamento del Meta (Hansen 2013)

A. Sitios detectados de quemas anuales (Giglio et al. 2003) para el departamento del Meta y B. Extensión de áreas deforestadas para el departamento del Meta (Hansen 2013)

# read data
##### ojo con def_meta
##### ojo con def_meta
##### ojo con def_meta
##### volver a version anteror en box que esta pegada con datos del SIGOT!!!!!
##### volver a version anteror en box que esta pegada con datos del SIGOT!!!!!
##### volver a version anteror en box que esta pegada con datos del SIGOT!!!!!

def_meta <- read.csv("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/municipi_table.csv")

# knitr::kable(head(def_meta), digits = 2, caption = "Deforestación por municipios")

def_meta2 <- read_excel ("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/db_final.xlsx")



# kable(head(def_meta))#, "html") %>%
  #kable_styling(bootstrap_options = c("striped","hover","condensed","responsive", full_width = F)) %>%
#scroll_box(width = "650px", height = "350px")

Fuegos

Los datos actividad de fuego provienen del sensor MODIS (Giglio et al. 2003). estos datos fueron obtenidos del servidor de la NASA https://data.nasa.gov/Earth-Science/MODIS-Active-Fire-Data/psbc-zh9i con un recorte para el departamento del Meta y agregados anualmente por municipio. Los datos estructurados en una matriz tipo panel y posteriormente fueron graficados por municipio para describir las tendencias de cambio, como un “proxy” antrópico en las sabanas naturales.

Variables explicativas

Se recopilaron y estructuraron datos de 1201 variables, para todos los municipios del Meta que fueron evaluados en cuanto a factores o causas de la transformación de ecosistemas. La base de datos se ensamblo de diversas fuentes como; DANE, DNP, Ministerio de Agricultura, IGAC, etc. Para ver una descripción detallada de las variables se puede consultar el siguiente documento [Documento descripción base de datos variables Meta (https://tnc.box.com/s/2zpwz3f5s86cimud21ywm8zv5xnb7suu)] y para ver la base de datos recopilada consulte [Tabla variables Meta (https://tnc.box.com/s/6rzr7sn3kbicnz1n5xfymo2ysothp7pt)].

Las 1201 variables fueron agrupadas en grandes grupos relacionados con el tipo de actividad agrícola, el tipo de uso y variables demográficas. Para agrupar las variables se realizó un análisis de correlación simple para seleccionar una o dos variables de cada grupo correlacionado.

Teniendo en cuenta el marco conceptual de (Geist and Lambin 2001) se seleccionaron variables que representaran los factores demográficos, económicos, políticos/institucionales y culturales/sociopolíticos. Una vez seleccionadas las variables en el marco conceptual de (Geist and Lambin 2001) se aplico el criterio de selección de variables que tuvieran el mismo rango de temporal de las variables de transformación de ecosistemas (2004-2010). El criterio rango temporal se aplico para poder realizar análisis que se aproximen más a la causalidad y no solamente una correlación espacial entre las variables explicativas y variables de tansformación de ecosistemas, ya que al detectar correlación espacial como temporal es más probable detectar una relación de causalidad. En el caso del factor tecnólogico no se encontraron variables que cumplieran con los criterios mínimos para su selección.

Diagrama de flujo de los análisis para explorar las causas de transformación de bosques y sabanas en el Meta y la caracterización socio-económica de los municipios

Diagrama de flujo de los análisis para explorar las causas de transformación de bosques y sabanas en el Meta y la caracterización socio-económica de los municipios

Las variables seleccionadas por los criterios del esquema conceptual de Geist 2001 y disponibles para los años 2000-2014 fueron:

  • Factor cultural y social: Resguardos Indigenas y cobertura educativa total.
  • Factor politico/institucional: Registro Único Nacional de Áreas Protegidas, acciones armadas y participación electoral.
  • Factor Demográfica: Población total
  • Factor Económico: Cultivos de coca, minería, costo de la tierra (valor catastral), actividad pecuaria y agrícola, vías, cultivos de palma y petroleo.

La variable de fuego fue integrada en el modelo para explicar la deforestación , con el objetivo de explorar si hay relación entre estos dos procesos de transformación de ecosistemas. Ya que en esta región las quemas usualmente estan asociadas a actividades para preparar las tierras para el uso agropecuario, pero al mismo tiempo las quemas pueden estar asociadas al proceso de deforestación (tumba de cobertura forestal y luego quema).

Posteriormente la base de datos fue analizada utilizando regresiones a pasos del modo (Backward), usando el criterio de informacion Akaike (AIC) como variable indicadora para dejar la variables fuera del modelo. El umbral sleccionado para el AIC fue de dos. De forma tal que el modelo inical considero todas las variables y al final nos quedamos con las covariables que explican la variacion estadistica en la deforestación y los fuegos. Ver figura 2

# read data set


datos_def <- def_meta

# datos_def$Deforest <- as.numeric(datos_def$Deforest)
# datos_def$Densidad <- as.numeric(datos_def$Densidad)

Selección de modelos

Usando estas variables filtradas, ajustamos varios modelos GLMM con los años y/o los municipios como “random-effect” para ver sus diferentes efectos. Los modelos más parsimoniosos se graficaron para simplificar.

Los GLMMs son una extensión de los modelos lineales generalizados (GLM), que incluyen efectos fijos y mezclados en la variable respuesta. Adicionalmente ofrecen una aproximación muy flexible para modelar datos un amplio rango de datos con estructuras complejas temporales y anidados (Harrison et al. 2018). Estos modelos siguen la forma general (en notación matricial):

\[y=X \beta + Z \gamma + \epsilon\]

Donde \(y\) es una columna de un vector N x 1, con la variable respuesta, en nuestro caso deforestación o fuego; X es una matriz N x p de p variables predictoras; \(\beta\) es una columna de un vector p x 1 de los efectos fijos de los coeficientes de la regresión (los “betas”); Z es una matriz de los efectos al azar multiplicado por \(\gamma\), el cual es un vector de los efectos al azar (el complemento a los efectos fijos en \(\beta\); y \(\epsilon\) es una columna o vector de los coeficientes residuales no explicados por el modelo.

\[\overbrace{\bf y}^{N\times 1}=\overbrace{\underbrace{\bf X}_{N\times p}\underbrace{\beta}_{p\times 1}}^{N\times 1}+\overbrace{\underbrace{\bf Z}_{N\times q}\underbrace{\gamma}_{q\times 1}}^{N\times 1}+\overbrace{\varepsilon}^{N\times 1}\]

Como ejemplo ilustrativo, si expandiéramos \(y\) por cada municipio en el modelo tendríamos:

\[y= \begin{vmatrix} deforest\\ 0.15 \\ 0.23 \\ \vdots \\ 0.45 \end{vmatrix}\]

y expandiendo X:

\[ X = \begin{vmatrix} Intercepto & avaluo & pozos & poblacion & ... & \\ 1 & 0.45 & 3 & 3500 & ...& \\ 0 & 0.55 & 5 & 1080 & ...\\ \vdots & \vdots & \vdots & \vdots & ...\\ 1 & 0.95 & 12 & 4200 & ...& \end{vmatrix}\]

Mientras que si expandimos \(\beta\)

\[ \beta= \begin{vmatrix} 0.22 \\ 0.35 \\ 0.33 \\ \vdots \\ 0.85 \end{vmatrix}\]

Dado Z es tan grande (29 municipios x 14 años) no describiremos su estructura. Mientras que \(\gamma\) es estimado asumiendo que se distribuye de forma normal como: \[{\gamma}\sim\mathcal{N}({ 0},{G})\] Donde G es una matriz simétrica de varianza-covarianza indexando la pendiente e intercepto.

library(ggplot2)
library(lme4)
library(nlme)
library(arm)
library(sjPlot)
library (readxl)# read excel data
library(MASS)
library(sjPlot)


def_meta2 <- as.data.frame(read_excel ("C:/Users/diego.lizcano/Box Sync/CodigoR/Cambio_Meta/data/db_final2.xlsx"))

def_meta2_avaluo <- as.numeric(def_meta2$avaluo)

def_meta2$educacion_total <- as.numeric(def_meta2$educacion_total)
def_meta2$avaluo <- as.numeric(def_meta2$avaluo)
def_meta2$eco_hidrocarb_mun_area_porc <- as.numeric(def_meta2$eco_hidrocarb_mun_area_porc)
def_meta2$eco_min_mun_area_porc <- as.numeric(def_meta2$eco_min_mun_area_porc)
def_meta2$pol_ar_mun_area_porc <- as.numeric(def_meta2$pol_ar_mun_area_porc)
def_meta2$pol_dcs_mun_area_porc <- as.numeric(def_meta2$pol_dcs_mun_area_porc)
def_meta2$pol_drmi_mun_area_porc <- as.numeric(def_meta2$pol_drmi_mun_area_porc)
def_meta2$pol_ley2_mun_area_porc <- as.numeric(def_meta2$pol_ley2_mun_area_porc)
def_meta2$pol_pnn_mun_area_porc <- as.numeric(def_meta2$pol_pnn_mun_area_porc)
def_meta2$pol_resgind_mun_area_porc <- as.numeric(def_meta2$pol_resgind_mun_area_porc)
def_meta2$pol_prn_mun_area_porc <- as.numeric(def_meta2$pol_prn_mun_area_porc)
def_meta2$pol_rfpn_mun_area_porc <- as.numeric(def_meta2$pol_rfpn_mun_area_porc)
def_meta2$pol_rnsc_mun_area_porc <- as.numeric(def_meta2$pol_rnsc_mun_area_porc)
def_meta2$fuego <- as.numeric(def_meta2$fuego)
def_meta2$pozos <- as.numeric(def_meta2$pozos)
def_meta2$loss <- as.numeric(def_meta2$loss)

def_meta2_poblacion <- as.numeric(def_meta2$poblacion)

# def_meta2[is.na(def_meta2)] <- 0
# def_meta2$avaluo <- def_meta2_avaluo
# def_meta2$poblacion <- def_meta2_poblacion

datos_scaled <- as.data.frame(scale(def_meta2[4:20], center = TRUE, scale = TRUE))
datos_scaled[is.na(datos_scaled)] <- 0
                     
datos_scaled$NOM_MUNICI <-  def_meta2$NOM_MUNICI               
datos_scaled$ano <-  def_meta2$ano

mod_glm1_forestloss <- glm (loss ~ ano + 
                  educacion_total +  avaluo + # no es cero es NA
                  eco_hidrocarb_mun_area_porc  +
                  eco_min_mun_area_porc +   
                  pol_drmi_mun_area_porc + 
                  pol_ley2_mun_area_porc +         
                  pol_pnn_mun_area_porc + 
                  pol_resgind_mun_area_porc +    
                  pol_prn_mun_area_porc +   
                  pol_rfpn_mun_area_porc +   
                  pol_rnsc_mun_area_porc + 
                  fuego + 
                  pozos + poblacion, # no es cero es NA
                   data=datos_scaled)

# Stepwise Regression

step_loss <- stepAIC(mod_glm1_forestloss, direction="backward", k = 2)
step_loss$anova # display results

mod_glm1_fire <- glm (fuego ~ ano + 
                              educacion_total +  avaluo + # no es cero es NA
                              eco_hidrocarb_mun_area_porc  +
                              eco_min_mun_area_porc +   
                              pol_drmi_mun_area_porc + 
                              pol_ley2_mun_area_porc +         
                              pol_pnn_mun_area_porc + 
                              pol_resgind_mun_area_porc +    
                              pol_prn_mun_area_porc +   
                              pol_rfpn_mun_area_porc +   
                              pol_rnsc_mun_area_porc + 
                              fuego + 
                              loss + poblacion, # no es cero es NA
                            data=datos_scaled)

step_fire <- stepAIC(mod_glm1_fire, direction="backward", k = 2)
step_fire$anova # display results


#### with Final Model:
# Model without respect to grouping
mod_glm2_forestloss <- glm (loss ~ avaluo + 
                              pol_rfpn_mun_area_porc + 
                              fuego + 
                              pozos + 
                              poblacion, 
                              data=datos_scaled)

summary(mod_glm2_forestloss)


# Model without respect to grouping
mod_glm2_fire <- glm (fuego ~ avaluo + 
                        eco_hidrocarb_mun_area_porc + 
                        pol_rfpn_mun_area_porc + 
                        pol_rnsc_mun_area_porc + 
                        loss + 
                        poblacion, 
                        data=datos_scaled)

summary(mod_glm2_fire)

#Observacción gráfica de la distribución de las variables 
# plot(step_loss)
# plot(step_fire)

# par(mfrow=c(5,1)) 
# hist(datos_scaled$avaluo)
# hist(datos_scaled$pol_rfpn_mun_area_porc)
# hist(datos_scaled$fuego)
# hist(datos_scaled$pozos)
# hist(datos_scaled$poblacion)
# dev.off()



#   MODELOS MIXTOS
#   Re-analizando los datos con modelos lineales mixtos: El paquete nlme




# 
# # Model with varying slope
# m2 <- lmer(salary ~ experience + (0 + experience|department), data=df)
# df$random.slope.preds <- predict(m2)
# 
# # Model with varying slope and intercept
# m3 <- lmer(salary ~ experience + (1 + experience|department), data=df)
# df$random.slope.int.preds <- predict(m3)
# 



# Model with varying intercept
library(Matrix)
library(nlme)


lmer1 <- lmer (loss ~ avaluo + 
               pol_rfpn_mun_area_porc + 
               fuego + 
               pozos + 
               poblacion +
              (1|NOM_MUNICI),  
               data=datos_scaled)


lmer2 <- lmer (loss ~ avaluo + 
               pol_rfpn_mun_area_porc + 
               fuego + 
               pozos + 
               poblacion +
               (1|ano), 
             data=datos_scaled)

# Model with varying slope
lmer3 <- lmer (loss ~ avaluo + 
                 pol_rfpn_mun_area_porc + 
                 fuego + 
                 pozos + 
                 poblacion +
                 (0+ano|NOM_MUNICI),#multiple uncorrelated random effects
               data=datos_scaled)



anova( lmer1, lmer2, lmer3)
summary(lme3)

plot_model(lmer3)
plot_model(lmer3, "re")
plot_model(lmer3, "slope")

sjp.lmer(lmer1,
         facet.grid = FALSE,
         sort.est = "sort.all",
         y.offset = .5, main="loss")


lmer4 <- lmer (fuego ~ avaluo + 
               eco_hidrocarb_mun_area_porc + 
               pol_rfpn_mun_area_porc + 
               pol_rnsc_mun_area_porc + 
               loss + 
               poblacion +
               (1|NOM_MUNICI), 
             data=datos_scaled )


# #  lme5 <- lme (fuego ~ avaluo + 
#                eco_hidrocarb_mun_area_porc + 
#                pol_rfpn_mun_area_porc + 
#                pol_rnsc_mun_area_porc + 
#                loss + 
#                poblacion +
#                (1|ano), 
#              data=datos_scaled )

# anova( lmer4, lmer5)

summary(lmer4)
# par(mfrow=c(1,2)) 
# plot(ranef(lme5), main="fuego - año")
# plot(ranef(lme4), main="fuego - municipio", vline=0)
plot_model(lmer4)
plot_model(lmer4, "re")
plot_model(lmer4, "slope")

sjp.lmer(lmer4,
         facet.grid = FALSE,
         sort.est = "sort.all",
         y.offset = .4)


# plot_model(lmer4, "eff", terms = "avaluo")

RESULTADOS

library(ggplot2)


# sjp.setTheme(legend.pos = "bottom")
             
sjp.scatter(datos_def$year,
            datos_def$loss, 
            datos_def$depto, 
            fit.line.grps = TRUE, 
            # fitmethod = "loess",
            # geom.colors = "Paired",
            fit.line = TRUE,
            show.ci = TRUE, 
            facet.grid = TRUE) #+ guides(nrow = 23)

Deforestación

Con el fin de explorar que factores y las variables que estan más correlacionadas con los patrones de deforestación durante el periodo 2000-2014 en los municipios del departamento del Meta. Las variables utilizadas para los modelos basado en el criterio de representar los factores propuestos por Lambin y Geist 2003 y que explicaran la variabilidad estadistica del fenomeno de la deforestación fueron: valor avaluo catastral, población total del municipio, frecuencia anual de fuegos, proporción del área del municipio con reserva forestal protectora y número de pozos.

La Figura 3 muestra los patrones espaciales (localización y extensión) para cada uno de los años que se realizó el análisis.

Áreas deforestación en municipios del departamento del Meta periodo 2000-2014 basado en Hansen et al. 2013

Áreas deforestación en municipios del departamento del Meta periodo 2000-2014 basado en Hansen et al. 2013

En la Figura 4 se muestran los coeficientes de los valores de deforestación del modelo de regresión para cada municipio. Un valor mayor a cero, en color azul, expresa una relación positiva en el modelo. En color rojo un valor de correlación negativo. Si la varianza no se cruza con la línea horizontal (cero) expresa una relación significativa. Se presentan en orden de la parte superior a la inferior los municipios de mayor a menor extensión deforestada durante el periodo 2000-2014. Los municipios como La Macarena, Puerto Rico la Uribe, Vista Hermosa, y Mesetas, presentan una reducción significativa de su áreas de bosques comparado con los municipios restantes del Meta (Figura 4). Los cuales contrastan con Puerto Gaitán, que presenta un incremento en la extensión de bosques.

Resultados de las asociaciones entre las variables explicativas respecto a la deforestación en el departamento del Meta. B. Resultado del mejor modelo que explica la relación entre las variables evaluadas y la deforestación en el departamento del Meta

Resultados de las asociaciones entre las variables explicativas respecto a la deforestación en el departamento del Meta. B. Resultado del mejor modelo que explica la relación entre las variables evaluadas y la deforestación en el departamento del Meta

En la Figura 5 se muestran los resultados obtenidos de los modelos donde se grafíca la relación entre cada una de las variables analizadas respecto al fenomeno de la deforestación donde se evidencia el tipo de asociación, la magnitud y la variabilidad. Las variables valor de avaluos catastrales por municipio, población total por municipio, proporción de áreas en reservas forestales protectoras naturales del municipio, número de pozos petroleros en el municipio y fuegos son buenos predictores del cambio en la cobertura forestal (Figura 5).

Los municipios con menores valores de avaluo catastral y menor presencia de pozos pretroléros presentan mayores tasas de deforestación. La variable de fuegos mostro una relación directamente proporcional con la deforestación. La varible de reservas forestales protectoras naturales muestra una relación positiva respecto a la deforestación sin embargo, esta relación presenta una variablidad alta. Y para el caso de población total mues tra una relación positiva sin embargo la magnitud de la relación con la deforestación es pequeña y no llega a ser significartiva.

Fuegos

Con el fin de explorar que factores y sus variables que estan más correlacionadas con los patrones de fuegos durante el periodo 2000-2014 en los municipios del departamento del Meta. Las variables utilizadas para los modelos basado en el criterio de representar los factores propuestos por Lambin y Geist 2003 y que explicaran la variabilidad estadistica de los fuegos fueron: valor avaluo catastral, población total del municipio, frecuencia anual de fuegos, proporción del área del municipio con bloques petroleros, proporción del área del municipio con reserva forestal protectora y proporción del área del municipio con reservas de la sociedad civil.

En la figura 6 se presentan los resultados del número de eventos de quemas por municipio, datos rescalados y acumulados para el periodo 2000-2014. Respecto a fuegos los municipios que muestran mayor número de eventos de fuego en el periodo analizado son: Puerto Gaitan, Mapiripán, San Martín y Puerto López. Y municipios como Puerto Rico y San Carlos de Guaroa presentan una baja frecuencía de eventos de fuegos.

Sitios de fuegos en municipios del departamento del Meta basado en Giglio et al. 2003

Sitios de fuegos en municipios del departamento del Meta basado en Giglio et al. 2003

Los resultados del modelo implementado resumidos en la figura 7 muestran una relación positiva entre los fuegos con las variables de avaluo catastral y reservas de la sociedad civíl. Sin embargo, la proporción de áreas de reservas de la sociedad civíl y reservas foretales naturales protectoras presentaron una alta variabilidad en su relación con los fuegos. Las variables de proporción de bloques petroleros y población total mostraron una relación inversamente proporcional con la variable de fuego.

Resultados de las asociaciones entre las variables explicativas respecto a quemas en el departamento del Meta. B. Resultado del mejor modelo que explica la relación entre las variables evaluadas y quemas en el departamento del Meta

Resultados de las asociaciones entre las variables explicativas respecto a quemas en el departamento del Meta. B. Resultado del mejor modelo que explica la relación entre las variables evaluadas y quemas en el departamento del Meta

Caracterización de variables socioeconómicas de los municipios del Meta

En el siguiente enlace pueden consultar las 1201 variables estudiadas para el análisis. La descripción completa de las variables y la base de datos se pueden encontrar en los siguientes enlaces (https://tnc.box.com/s/2zpwz3f5s86cimud21ywm8zv5xnb7suu) y (https://tnc.box.com/s/6rzr7sn3kbicnz1n5xfymo2ysothp7pt).

Valores de avaluao catastral por municipios del Meta para el año 2016

Valores de avaluao catastral por municipios del Meta para el año 2016

CONCLUSIONES

Los municipios con mayor deforestación en el Meta están más asociados a bajos valores de costo de la tierra (avalúo) y menor presencia de actividades económicas como el petróleo (número de pozos). Y sugiere que la deforestación asociada a áreas donde el valor de la tierra es bajo y donde posiblemente se dan actividades informales de baja escala. El aumento de extensiones de bosques en municipios como Puerto Gaitán, San Martin y San Luis de Cubarral puede responder a la implementación de sistemas forestales (por ejemplo, caucho y palma entre otros) que no corresponde al incremento en la extensión bosques naturales.

En relación a los municipios con mayor frecuencia de eventos de fuegos, se encontró están más asociados a tierras con mayor valor económico (avalúo) y municipios con menor población total. Estos patrones pueden estar asociados a actividades económicas de ganadería y agricultura donde las prácticas de quemas suelen estar asociadas a estas actividades. Pese a que los datos distribuidos en por Giglio y colaboradores en teoría remueven las detecciones de fuego por teas asociadas a la producción petrolera, se recomienda revisar si los patrones de frecuencia y posición de las quemas excluyen esta actividad ya que en el caso de que no sea así esta actividad no necesariamente muestra una progresión en la transformación de ecosistemas naturales.

Los resultados muestran que los patrones de quemas y deforestación son fenómenos que se encuentran relacionados, y soporta los patrones encontrados por múltiples estudios donde estas dos actividades están relacionadas en el proceso de preparar los suelos para el desarrollo de actividades pecuarias y agrícolas. En este sentido las quemas pueden ser una actividad asociada a la transición de sabanas y bosques hacia otras coberturas y usos de la tierra.

En cuanto a las variables utilizadas en los análisis, son pocas las variables disponibles que cumplieran con los criterios seleccionados para realizar análisis robustos que se aproximan a una relación de causalidad, por lo que los resultados acá presentados son una exploración en la asociación de variables socioeconómicas con quemas y deforestación en el departamento del Meta bajo este enfoque. El entendimiento de las causas próximas y subyacentes requieren de análisis a diferentes escalas espaciales y temporales, donde los datos que se integren cumplan con criterios de selección que permitan encontrar patrones asociados a la pregunta y es importante generar/seleccionar variables bajo criterios que estén acorde a los supuestos de los análisis.
Como recomendaciones de pasos a seguir se sugiere explorar métodos alternativos y complementarios para analizar los cambios de uso y cobertura de la tierra y factores asociados a estos, a múltiples escalas temáticas, espaciales y temporales. El conocimiento de las causas de la transformación de ecosistemas es un proceso dinámico que integra múltiples disciplinas y fuentes de conocimiento, en ese sentido el presente trabajo sintetiza una aproximación estadística que contribuye al entendimiento de las dinámicas de cambio del Meta. Se sugiere que para futuros ejercicios se debe partir de las preguntas que se quieren responder con el estudio de causas de transformación para definir los métodos e información a utilizar.

Los resultados de caracterización socioeconómica de los municipios del Meta proveen información complementaria a este análisis de causas de transformación para que usuarios y múltiples análisis más allá de este análisis exploren e integren información sistematizada y estructurada a futuros análisis y preguntas.

Anexos

Agrupando variables correlacionadas

A partir de una tabla con 1201 variables que fueron recopiladas, se exploraron correlaciones entre grupos variables, para eliminar variables correlacionadas.

Variables agropecuarias

# Correlation matrix from info_meta
# with mpg, cyl, and disp as rows 
# and hp, drat, and wt as columns 
x1 <- info_meta[151:200]
y1 <- info_meta[151:200]
corMat1 <- cor(x1, y1)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat1, order = "FPC",type = "lower", tl.pos = "lower", tl.cex = 0.5)

Variables agropecuarias-tendencia

x2 <- info_meta[201:230]
y2 <- info_meta[201:230]
corMat2 <- cor(x2, y2)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat2, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuaria-uso

x3 <- info_meta[231:255]
y3 <- info_meta[231:255]
corMat3 <- cor(x3, y3)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat3, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-vivienda

x4 <- info_meta[262:285]
y4 <- info_meta[262:285]
corMat4 <- cor(x4, y4)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat4, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-maquinaria

x5 <- info_meta[292:333]
y5 <- info_meta[292:333]
corMat5 <- cor(x5, y5)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat5, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-infraestructura

x6 <- info_meta[334:360]
y6 <- info_meta[334:360]
corMat6 <- cor(x6, y6)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat6, order = "original",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-asistencia tecnica

x7 <- info_meta[361:396]
y7 <- info_meta[361:396]
corMat7 <- cor(x7, y7)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat7, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-agua

x8 <- info_meta[397:465]
y8 <- info_meta[397:465]
corMat8 <- cor(x8, y8)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat8, order = "original",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-suelo

x9 <- info_meta[466:495]
y9 <- info_meta[466:495]
corMat9 <- cor(x9, y9)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat9, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-ganado

x10 <- info_meta[541:594]
y10 <- info_meta[541:594]
corMat10 <- cor(x10, y10)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat10, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Variables agropecuarias-porcinos

x10 <- info_meta[595:630]
y10 <- info_meta[595:630]
corMat11 <- cor(x10, y10)

# corrplot(corMat, method = "color", type="lower")#,
         # col = colorRampPalette(c("blue", "white", "red"))(100))

corrplot(corMat11, order = "FPC",type = "lower", tl.pos = "full", tl.cex = 0.5)

Session Info

Detalles y paquetes de la sesión en R

print(sessionInfo(), locale = FALSE)
## R version 3.4.0 (2017-04-21)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
## 
## Matrix products: default
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_0.8.0    knitr_1.20          corrplot_0.84      
##  [4] printr_0.1          gfcanalysis_1.4     sf_0.6-3           
##  [7] readxl_1.1.0        mapview_2.4.0       forcats_0.3.0      
## [10] stringr_1.3.1       dplyr_0.7.4         purrr_0.2.4        
## [13] readr_1.1.1         tidyr_0.8.1         tibble_1.4.2       
## [16] tidyverse_1.2.1     ggmap_2.6.1         ggplot2_2.2.1      
## [19] rgdal_1.2-20        rgeos_0.3-26        maptools_0.9-2     
## [22] rasterVis_0.44      latticeExtra_0.6-28 RColorBrewer_1.1-2 
## [25] lattice_0.20-35     raster_2.6-7        sp_1.2-7           
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.3-2   rjson_0.2.18       class_7.3-14      
##  [4] gdalUtils_2.0.1.14 leaflet_2.0.0      rprojroot_1.3-2   
##  [7] satellite_1.0.1    base64enc_0.1-3    rstudioapi_0.7    
## [10] hexbin_1.27.2      lubridate_1.7.4    xml2_1.2.0        
## [13] codetools_0.2-15   R.methodsS3_1.7.1  mnormt_1.5-5      
## [16] jsonlite_1.5       broom_0.4.4        png_0.1-7         
## [19] R.oo_1.22.0        RPostgreSQL_0.6-2  shiny_1.1.0       
## [22] mapproj_1.2.6      compiler_3.4.0     httr_1.3.1        
## [25] backports_1.1.2    assertthat_0.2.0   lazyeval_0.2.1    
## [28] cli_1.0.0          later_0.7.2        htmltools_0.3.6   
## [31] tools_3.4.0        bindrcpp_0.2.2     gtable_0.2.0      
## [34] glue_1.2.0         reshape2_1.4.3     maps_3.3.0        
## [37] Rcpp_0.12.16       cellranger_1.1.0   nlme_3.1-137      
## [40] udunits2_0.13      iterators_1.0.9    crosstalk_1.0.0   
## [43] psych_1.8.4        proto_1.0.0        rvest_0.3.2       
## [46] mime_0.5           zoo_1.8-1          scales_0.5.0.9000 
## [49] hms_0.4.2          promises_1.0.1     parallel_3.4.0    
## [52] animation_2.5      yaml_2.1.19        geosphere_1.5-7   
## [55] stringi_1.2.2      foreach_1.4.4      e1071_1.6-8       
## [58] spData_0.2.8.3     RgoogleMaps_1.4.1  rlang_0.2.0       
## [61] pkgconfig_2.0.1    bitops_1.0-6       evaluate_0.10.1   
## [64] bindr_0.1.1        htmlwidgets_1.2    plyr_1.8.4        
## [67] magrittr_1.5       R6_2.2.2           DBI_1.0.0         
## [70] pillar_1.2.2       haven_1.1.1        foreign_0.8-67    
## [73] units_0.5-1        RCurl_1.95-4.10    modelr_0.1.2      
## [76] crayon_1.3.4       rmarkdown_1.9      jpeg_0.1-8        
## [79] grid_3.4.0         digest_0.6.15      classInt_0.2-3    
## [82] webshot_0.5.0      xtable_1.8-2       httpuv_1.4.3      
## [85] R.utils_2.6.0      stats4_3.4.0       munsell_0.4.3     
## [88] viridisLite_0.3.0

Literatura Citada

Acevedo-Charry, O., A. Pinto-Gómez, and J. O. Rangel-Ch. 2014. Las aves de la orinoquia colombiana: Una revisión de sus registros. Colombia Diversidad Biótica XIV. La región de la Orinoquía de Colombia:691–750.

Armenteras, D., E. Cabrera, N. Rodríguez, and J. Retana. 2013a. National and regional determinants of tropical deforestation in Colombia. Regional Environmental Change.

Armenteras, D., N. Rodríguez, and J. Retana. 2013b. Landscape Dynamics in Northwestern Amazonia: An Assessment of Pastures, Fire and Illicit Crops as Drivers of Tropical Deforestation. PLoS ONE 8:e54310.

Armenteras, D., N. Rodríguez, and J. Retana. 2013c. Landscape dynamics in northwestern Amazonia: an assessment of pastures, fire and illicit crops as drivers of tropical deforestation. PloS one 8:e54310.

Armenteras, D., N. Rodríguez, J. Retana, and M. Morales. 2011. Understanding deforestation in montane and lowland forests of the Colombian Andes. Regional Environmental Change 11:693–705.

Boron, V., J. Tzanopoulos, J. Gallo, J. Barragan, L. Jaimes-Rodriguez, G. Schaller, and E. Payán. 2016. Jaguar Densities across Human-Dominated Landscapes in Colombia: The Contribution of Unprotected Areas to Long Term Conservation. PLOS ONE 11:e0153973.

Brooks, T. M., R. A. Mittermeier, C. G. Mittermeier, G. A. B. da Fonseca, A. B. Rylands, W. R. Konstant, P. Flick, J. Pilgrim, S. Oldfield, G. Magin, and C. Hilton-Taylor. 2002. Habitat loss and extinction in the hotspots of biodiversity. Conservation Biology 16:909–923.

Etter, A., C. McAlpine, and H. Possingham. 2008. Historical Patterns and Drivers of Landscape Change in Colombia Since 1500: A Regionalized Spatial Approach. Annals of the Association of American Geographers 98:2–23.

Foley, J. A., R. Defries, G. P. Asner, C. Barford, G. Bonan, S. R. Carpenter, F. S. Chapin, M. T. Coe, G. C. Daily, H. K. Gibbs, J. H. Helkowski, T. Holloway, E. A. Howard, C. J. Kucharik, C. Monfreda, J. A. Patz, I. C. Prentice, N. Ramankutty, and P. K. Snyder. 2005. Global consequences of land use. Science (New York, N.Y.) 309:570–4.

Furley, P. A. 1999. The nature and diversity of neotropical savanna vegetation with particular reference to the brazilian cerrados. Global Ecology and Biogeography 8:223–241.

Geist, H. J., and E. F. Lambin. 2001. What drives tropical deforestation. LUCC Report series 4:116.

Giglio, L., J. Descloitres, C. O. Justice, and Y. J. Kaufman. 2003. An Enhanced Contextual Fire Detection Algorithm for MODIS. Remote Sensing of Environment 87:273–282.

Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice, and J. R. G. Townshend. 2013. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 342:850–853.

Harrison, X. A., L. Donaldson, M. E. Correa-Cano, J. Evans, D. N. Fisher, C. E. Goodwin, B. S. Robinson, D. J. Hodgson, and R. Inger. 2018. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6:e4794.

IDEAM, I. d. 2010. Estudio nacional del agua 2010. Bogotá DC.

Jackson, C. H. 2011. Multi-State Models for Panel Data: The msm Package for R. J Stat Softw 38:1–28.

Jha, S., and K. S. Bawa. 2006. Population growth, human development, and deforestation in biodiversity hotspots. Conservation Biology 20:906–912.

Lambin, E. F., H. J. Geist, and E. Lepers. 2003. Dynamics of land-use and land-cover change in tropical regions. Annual Review of Environment and Resources 28:205–241.

Maldonado-Ocampo, J. A., A. Urbano-Bonilla, J. V. Preciado, and J. D. Bogotá-Gregory. 2013. Peces de la cuenca del río pauto, orinoquia colombiana. Biota colombiana 14.

Romero, M., J. Maldonado, S. Usama, A. M. Umaña, J. Murillo, S. Restrepo, M. Alvarez, M. T. Palacios, S. Valbuena, S. L. Mejía, and others. 2009. Informe sobre el estado de la biodiversidad en colombia 2007-2008: Piedemonte orinoquense, sabanas y bosques asociados al norte del río guaviare. Romero, MH et al. 2009. Informe sobre el estado de la biodiversidad en Colombia 2007-2008: Piedemonte orinoquense, sabanas y bosques asociados al norte del río Guaviare Informe sobre el estado de la biodiversidad en Colombia 2007-2008: Piedemonte orinoquense, sabanas y bosques asociados al norte del río Guaviare. Bogotá: Instituto de Investigación de Recursos Biológicos Alexander von Humboldt. 154 p.

Romero-Ruiz, M., S. Flantua, K. Tansey, and J. Berrio. 2012. Landscape transformations in savannas of northern South America: Land use/cover changes since 1987 in the Llanos Orientales of Colombia. Applied Geography 32:766–776.

Sánchez-Cuervo, A. M., T. M. Aide, M. L. Clark, and A. Etter. 2012. Land cover change in Colombia: surprising forest recovery trends between 2001 and 2010. PloS one 7:e43943.

Vitousek, P. M., H. A. Mooney, J. Lubchenco, and J. M. Melillo. 1997. Human domination of earth’s ecosystems. Science 277:494–499.