Documento1_climaJBM.Rnw

\documentclass[spanish,12pt]{article}
\usepackage[utf8]{inputenc} 
\usepackage[T1]{fontenc}
\usepackage[spanish]{babel}
\usepackage{csquotes}
\usepackage{lscape}

\usepackage[shiftmargins]{vmargin}
\setpapersize{USletter}
%%\setmarginsrb{left}{top}{right}{bottom}{headhgt}{headsep}{foothgt}{footskip}
\setmarginsrb{3.2cm}{2cm}{1.8cm}{2cm}{0.5cm}{0.5cm}{0.5cm}{0.5cm}

\usepackage[Glenn]{fncychap}
\usepackage{longtable}
%% cambiamos a biblatex
%%\usepackage[backend=bibtex,natbib=true,style=authoryear]{biblatex}
%\addbibresource{/home/jferrer/CEBA/lib/BibTEX/JBM.bib}

\usepackage[authoryear,sort]{natbib}



\usepackage[pagebackref=false,colorlinks=true,citecolor=black,linkcolor=black,filecolor=black,urlcolor=black]{hyperref}

\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\usepackage{titling}

\pretitle{\begin{center} 
\large \colorbox{black}{\textcolor{white}{Laboratorio de Ecología Espacial}} \par \HRule \\[0.4cm] \LARGE }
\posttitle{\par \HRule \\[1.5cm]
%\includegraphics[width=8cm]{/Users/jferrer/NeoMapas/img/logo_NeoMapas.jpg}
\end{center}\vskip 0.5em}

\preauthor{\large\begin{center}}
\postauthor{\\
Laboratorio de Ecología Espacial, Centro de Estudios Botánicos y Agroforestales\\Instituto Venezolano de Investigaciones Científicas
\end{center}
\par}

\predate{\begin{center}\normalsize Borrador\\\Sexpr{cdg.doc}\\}
\date{Versión de \today}
\postdate{\par\end{center}}


\title{Estacionalidad y tendencias en el clima y la vegetación de un fragmento semi-urbano de bosque seco en el municipio San Francisco, estado Zulia}
\author{JR Ferrer-Paris}

\begin{document}

\maketitle

\bibliographystyle{/home/jferrer/NeoMapas/lib/custom-bib/tesis}


\section*{A quien pueda interesar}

<<echo=false, eval=true, results=hide>>=
require(chron)
require(raster)
require(RColorBrewer)

paquetes <- (.packages())
paquetes <- paquetes[!(paquetes %in% c("stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "deldir", "DBI", "RMySQL"))]

luq <- function(x,contar.NA=FALSE) {
	if (contar.NA==F) {
	x <- x[!is.na(x)]
	}
 length(unique(x))
 }

@

Este documento es generado utilizando las funciones de \texttt{Sweave} desde una sesión de \emph{R} \citep{pqt::base}, por tanto todas las tablas y figuras se generan y actualizan automáticamente a partir de los datos suministrados. Para acceso al código fuente en \emph{R} y los archivos de datos contacte al primer autor. Dentro de \emph{R} utilizamos los paquetes 
<<citas paquetes, results=tex, echo=false, eval=true>>=
cat(paste("\\emph{",paquetes,"} \\citep{pqt::",paquetes,"}",sep="",collapse="; "))
@
.

\section{Introducción}

Los bosques secos de la cuenca del Lago de Maracaibo han sufrido una alta tasa de deforestación, fragmentación y transformación y se estima que sólo se mantiene el 12\% de la cobertura original \citep{Portillo2012}. Uno de los fragmentos remanentes más emblemáticos está ubicado en el sector conocido como La Sibucara o Caujarito, a unos 20km al sur de la ciudad de Maracaibo entre los terrenos del Jardín Botánico de Maracaibo, la Base Aérea Rafael Urdaneta y el Aeropuerto Internacional “La Chinita”. La situación de este fragmento es vulnerable pues tiene una extensión reducida, y está rodeado por una matriz de actividades intensivas (viviendas, infraestructura, ganadería, etc.), pero a la vez tiene un valor estratégico por su ubicación dentro de la Zona Protectora de Maracaibo, entre los municipios Maracaibo y San Francisco, y con un gran potencial de uso recreacional, educativo y científico \citep{reporteZPM}. 

No es posible establecer la extensión original de la vegetación boscosa natural, pero existe evidencia de que esta área estuvo sometida a explotación por pastoreo y otros usos durante muchos años \citep{Martinez1968,Zambrano77}. En la década de 1950 se inicia la construcción del Aeropuerto Internacional La Chinita, el cual sería inaugurado oficialmente en marzo de 1969. A partir de 1975 se inician las operaciones en la Base Aérea Rafael Urdaneta \citep{HernandezParra99}.

Los reglamentos aeronáuticos establece una serie de restricciones en las actividades y edificaciones que se pueden realizar en sus inmediaciones, los cuales son reafirmados en el reglamento de uso de la Zona Protectora del Área Metropolitana de la ciudad de Maracaibo \citep{ReglamentoZPM}. Esta situación ha limitado la intervención humana dentro de este fragmento en los últimos 40 años, permitiendo el desarrollo y regeneración de la vegetación natural. Sin embargo, la cara norte de este fragmento constituye una excepción notable. 

Durante los años 1979 a 1983 se produjo una intervención importante por la creación del Jardín Botánico de Maracaibo (JBM), inaugurado oficialmente en 1982, que alteraría el paisaje natural e introduciría una serie de especies vegetales exóticas. Durante gran parte de los años 80 estuvo considerado entre los 8 principales jardínes botánicos de Venezuela, entre los que destacaba por su programa pionero de formación en horticultura, pero la continuidad de sus operaciones sufrió crónicamente por la falta de apoyo institucional y la declinación progresiva de personal técnico y profesional \citep{Huber1998}. Para 1992 habían cesado casi por completo las actividades básicas de mantenimiento y cuidado de las áreas verdes \citep{HernandezParra99}. A partir de esta fecha se observaron procesos de transformación natural de la vegetación en las áreas abandondas, con intervenciones puntuales por actividades de baja intensidad realizadas por vecinos y residentes (extracción forestal, pastoreo y cultivos focalizados, \citetext{Chourio, N., com. pers.}). 

Los esfuerzos de la Fundación Jardín Botánico de Maracaibo por recuperar el área con el apoyo de instituciones públicas y privadas se concretaron entre 2012 y 2013 con un Plan de Rehabilitación Integral del JBM. En estos años se realizaron obras de intervención y limpieza para la habilitación básica y mantenimiento de 20 hectáreas de jardín y la rehabilitación de la infraestructura de servicios y algunas estructuras recreativas que permitieron la apertura parcial acompañada de diversas labores de difusión y promoción \citep{InformeFJBM}.

El contraste entre la vegetación natural de las inmediaciones del aeropuerto y la vegetación manejada, abandonada y posteriormente recuperada dentro del JBM puede aportar datos útiles para futuros esfuerzos de recuperación y restauración en ecosistemas secos de la cuenca del lago de Maracaibo. Entre 1973 y 1974 se realizó un estudio de la vegetación de la zona y se desarrolló una florula local para caracterizar el bosque antes de la intervención que significó el desarrollo del JBM \citep{Zambrano77}. Lamentablemente no se han encontrado documentos o estudios posteriores que permitan evaluar el impacto de la introducción de especies y los cambios en la composición florística de las áreas intervenidas durante el periodo inicial de funcionamiento del JBM, y aún no se ha realizado una revisión y evaluación de las muestras recuperadas del herbario del JBM, actualmente depositadas en el Herbario de la Facultad de Agronomía, HERZU \citep{InformeFJBM}. Actualmente diversas instituciones, entre ellas La Universidad del Zulia (LUZ) y el Instituto Venezolano de Investigaciones Científicas (IVIC), están promoviendo estudios locales sobre la cobertura vegetal, y la fauna de vertebrados de la zona, entre otros \citep{Fernandez2011,Gonzalez2012,Larreal2012,Quintero2013,Ortega2013,Cardozo2014}.

El objetivo del presente estudio es caracterizar el clima y la fenología de la vegetación en este importante fragmento boscoso a partir de datos de estaciones climáticas y sensores remotos, a fin de sentar una linea base para estudios futuros. Los datos recabados permiten señalar tendencias temporales previas al año 2012, y que pueden ayudar a entender las transformación que ha sufrido la zona en el pasado reciente y que podría sufrir en el futuro cercano. 

\section{Métodos}

\subsection{Área de estudio}

%\citet{Zambrano77} clasifican el área como bosque muy seco tropical.  

%En el mapa de Vegetación de Huber y Alarcón (1988), esta área entra en la región B (Llanuras), subregión B1 (Depresión de Maracaibo), y pueden observarse formaciones de arbustales xerófilos espinosos ('cardonales y espinares', comunidades bajas entre 3-8m de altura, con gran proporción de especies armadas)  y bosques xerófilos bajos, deciduos (hasta 12-15m de altura, con dos estratos arbóreos deciduos y un sotobosque bien definido y denso). Esta región pertenece a la Provincia fitogeográfica del Caribe Meridional, pero aparentemente no presenta ninguna especie característica del denominado Distrito fitogeográfico Zuliano.



<<mapoteca,echo=false>>=
mptc <- "~/CEBA/lib/mapoteca/"
mapaDB <- "Maracaibo"
mapaDB <- "ZPM"

if (!exists("bosque")) {
  LC <- raster(sprintf("%s/%2$s/MCD12Q1.A2005001.%2$s_LC1.tif",mptc,mapaDB))
  bosque <- raster(sprintf("%s/%s/GFC2013/GFC2013.%2$s.treecover2000.tif",mptc,mapaDB))
  perdida <- raster(sprintf("%s/%s/GFC2013/GFC2013.%2$s.loss.tif",mptc,mapaDB))
  ganancia <- raster(sprintf("%s/%s/GFC2013/GFC2013.%2$s.gain.tif",mptc,mapaDB))
  frags <- clump(bosque>40)
  sfrags <- clump(ganancia)

  tt <- rev(sort(table(values(frags))))
  values(frags)[values(frags) %in% as.numeric(names(tt)[tt<200])] <- 9999
  ##table(values(sfrags))
  ##plot(sfrags %in% 446)
  tt <- rev(sort(table(values(frags))))
  
  slc <- c(426, 441, 442) ## para Maracaibo
  slc <- c(602, 617, 1504,449,98,277,1288) ## para ZPM
  ##plot(frags %in% slc[1:2])
  ##plot(frags %in% slc[c(4,6)])

##  r0 <- (frags %in% slc[1:2]) + sfrags %in% 446 ## para Maracaibo
    r0 <- (frags %in% slc[1:2]) + sfrags %in% 669 ## para ZPM
  JBM <- rasterToPolygons(r0, fun=function(x) {x %in% 1},dissolve=T)
  lindero.JBM <- shapefile("~/CEBA/data/JardinBotanicoMaracaibo/JBM.shp")
  rJBM <- rasterize(lindero.JBM,r0)
}

##area bosque fragmento
##sum(values(r0*area(r0)),na.rm=T)

## area del jardin
##sum(values(rJBM*area(r0)),na.rm=T)

## area bosque dentro del jardin
##sum(values(rJBM*r0*area(r0)),na.rm=T)

@ 


<<municipios,echo=false,results=hide>>=
vptc <- "~/NeoMapas/lib/gisdata/"                              
if (!exists("rADM2")) {
  
##  municipio <- shapefile(sprintf("%s/%s",mptc,"VEN_adm2.shp"))
##  subset(municipio@data,NAME_2 %in% c("Maracaibo","San Francisco") & NAME_1 %in% "Zulia")
##  rADM2 <- rasterize(municipio,bosque,field="ID_2")
  
  
  municipio <- shapefile(sprintf("%s/%s",vptc,
                                 "Municipios Vzla_CGeog_WGS84_region.shp"))
  vialidad <- shapefile(sprintf("%s/%s",mptc,"VEN_roads.shp"))
  
  subset(municipio@data,iconv(MUNICIPIO,"latin1","utf8") %in% toupper(c("Maracaibo","San Francisco","La Cañada de Urdaneta","Mara")) & ESTADO %in% "ZULIA")

  rADM2 <- rasterize(municipio,bosque,field="CODIGO")

}
  xarea <- median(values(area(rADM2)))

@ 


<<echo=false,fig=true>>=
if (!exists("rZPM")) {
  ZPM <- read.csv(file="~/CEBA/data/JardinBotanicoMaracaibo/ZonaProtectoraCoords.csv",row.names=1)
  ZPM <- rbind(ZPM,head(ZPM,1))
  coordinates(ZPM) <- 1:2
  proj4string(ZPM) <- "+proj=utm +zone=19n +datum=WGS84"
  
  ZPM.ll <- spTransform(ZPM,bosque@crs)
  ZPM.ll <- SpatialPolygons(list(Polygons(list(Polygon(ZPM.ll)),ID=1)),proj4string=ZPM.ll@proj4string)
  
  rZPM <- rasterize(ZPM.ll,bosque)
}

plot(ZPM.ll,col="grey88",lty=3)
##plot(vialidad,add=T,col="grey",lwd=3)
##plot(vialidad,add=T,col=1,lty=2,lwd=1)
plot(municipio,lty=1,lwd=2,add=T)
plot(JBM,col="darkgreen",add=T)

g <- floor(abs(JBM@bbox))
g <- g*sign(JBM@bbox)
m <-  (abs(JBM@bbox) - floor(abs(JBM@bbox))) * 60
##floor(m)
s <-  (abs(m) - floor(abs(m))) * 60

xys <- sprintf("%s°%s'%0.2f''",g,floor(m),s)

@ 

El área de estudio está ubicada entre \Sexpr{xys[2]} y \Sexpr{xys[4]} de latitud norte y \Sexpr{xys[1]} y \Sexpr{xys[3]} de longitud oeste, en la Depresión de Maracaibo \citep[región B, subregión B1, del mapa de Vegetación de][]{HuberAlarcon1988}. La vegetación consiste en arbustales xerófilos espinosos y bosques xerófilos bajos, deciduos, pertenecientes a la provincia fitogeográfica del Caribe Meridional, dentro del denominado distrito fitogeográfico zuliano. 

%Para el año 2000 se estima que la cobertura boscosa en la ZPM es de ... km2 de los cuales 33,8 \% se encuentran en el municipio San Francisco, y el parche de bosque representa \% del total. 

<<Area,echo=false>>=
calc.area <- table(values(bosque)>40,values(rADM2),values(rZPM))
area.2000 <- sum(values(frags %in% slc[1:2]))*xarea
area.2012 <- sum(values(r0))*xarea

@ 
<<Area2,eval=false,echo=false>>=

table(values(bosque)>0,values(rADM2))
table(values(rZPM),values(rADM2))

area.2000 <- sum(values(rADM2) %in% c(275,311) & values(frags %in% slc[1:2]))*xarea
area.2012 <- sum(values(r0))*xarea

table(values(rADM2) %in% c(275,311),values(frags %in% slc[1:2]))*xarea
table(values(rADM2) %in% c(275,311),values(frags>0))*xarea

table(values(rADM2),values(frags))*xarea
round(table(values(rADM2),values(frags))*xarea,3)
##plot(bosque)
##plot(municipio,add=T)

table(values(frags),values(perdida))
table(values(sfrags %in% 446))*xarea
@ 

Estimados recientes de cobertura boscosa basada en imágenes de satélite LandSat sugieren que para el año 2000 el área de estudio contaba con dos fragmentos que sumaban \Sexpr{round(area.2000,3)} km2 o \Sexpr{round(area.2000*100,1)} hectáreas con porcentaje de cobertura arborea mayor a 40\% \citep{Hansen2013}.%, los cuales representan los dos fragmentos de cobertura boscosa de mayor tamaño en la Zona Protectora de Maracaibo.
En el periodo 2000 a 2012 se produjo un aumento notable de la cobertura boscosa en el área del JBM, y se incrementó el área total de los fragmentos en un \Sexpr{round((area.2012-area.2000)*100/area.2000,1)}\%. %en el área d, y el est  en el área total de los fragmentos, y pérdidas marginales de menos del 1%.
Trazamos un polígono con un área total de \Sexpr{round(area.2012*100,1)} \emph{ha} alrededor de estos fragmentos para delimitar el área de estudio

%En el 2000 el area estimada es de . Para el 2012, el área incrementó hasta \Sexpr{round(area.2012,3)} km2 o \Sexpr{round(area.2012*100,1)} ha.
%, y aproximadamente el 33,6% de la cobertura boscosa total, 

Falta colocar parroquia y municipio en la ubicación, linderos, etc...

\subsection{Datos de estación meteorológica}

La estación meteorológica 804070 (también referida como 80407, o con el serial nacional 1015) está ubicada en el Aeropuerto Internacional La Chinita en Maracaibo (N10°34' W71°44', 65 msnm), dentro del área de estudio. La estación cuenta con un registro histórico desde 1959 hasta el 2013, pero los datos previos a 1973 provienen de una ubicación diferente (antiguo Aeropuerto Grano de Oro). 

Para el presente análisis considero la serie de tiempo de mediciones diarias desde el primero de junio de 1973 hasta 31 de diciembre de 2013, descargada a través del enlace \url{http://www.tutiempo.net/clima/Maracaibo-La_Chinita/804070.htm}, y verificada con los promedios y totales anuales publicados para el periodo 1961 a 1990 \citep{NormClim1993}. La serie presenta valores vacíos en las mediciones, de temperatura y precipitación para 17.8\% y 21.6\% de los días, respectivamente, especialmente en la segunda mitad de la década de los 80. 

Para el análisis de la estacionalidad utilizo el los datos del periodo 2000 a 2011, agregados en periodos de 16 días comparables con los datos de la serie de tiempo de mediciones realizadas con sensores remotos. 

<<bajarDatosClimaticos,eval=false,echo=false>>=
setwd(" ~/CEBA/data/JardinBotanicoMaracaibo/climaviejo")
for (yy in 1959:2013) {
  for (mm in 1:12) {
    system(sprintf("wget --continue 'http://www.tutiempo.net/clima/Maracaibo-La_Chinita/%02d-%04d/804070.htm' --output-document=/home/jferrer/CEBA/data/JardinBotanicoMaracaibo/climaviejo/E804070_%04d_%02d.htlm",mm,yy,yy,mm))
  }
}

@ 




<<DatosClimaticos,echo=false>>=
if (file.exists("~/CEBA/Rdata/JBMts.rda")) {
  load(file="~/CEBA/Rdata/JBMts.rda")
}
if (!exists("dts.clm")) {
  dts.clm <- data.frame()
  for (year in 1959:2013) {
    for (month in 1:12) {
      if (file.exists(sprintf("~/CEBA/data/JardinBotanicoMaracaibo/climaviejo/E804070_%i_%02i.htlm",year,month))) {
        system(sprintf("html2text ~/CEBA/data/JardinBotanicoMaracaibo/climaviejo/E804070_%i_%02i.htlm > prueba.txt",year,month))
      
        ##j <- as.numeric(system("grep -n 'Valores medios climáticos' prueba.txt | cut -d: -f1",intern=T))+2
        j <- as.numeric(system("grep -n 'Valores históricos' prueba.txt | cut -d: -f1",intern=T))+1
        k <- as.numeric(system("grep -n 'Medias y totales mensuales' prueba.txt | cut -d: -f1",intern=T))-1
        if (length(j)>0) {
          system(sprintf("sed -n %s,%sp prueba.txt > tabla",j,k))
          ##k=`grep -n "Medias y totales mensuales" prueba.txt| cut -d: -f1`
          ##let "j += 2"
          ##let "k -= 1"
          ##sed -n $j,$kp prueba.txt > tabla
          ##sed -n $j,+31p prueba.txt > tabla
          
          
          tt <- read.table("tabla",sep="",as.is=T)
          dts.clm <- rbind(dts.clm,data.frame(year,month,day=1:nrow(tt),
                                        tmean=as.numeric(tt$DíaT),
                                        tmin=as.numeric(tt$Tm),
                                        tmax=as.numeric(tt$TM),
                                        H=as.numeric(tt$H),
                                        PP=as.numeric(tt$PP),
                                        VV=as.numeric(tt$VV)))
        }
      }
    }
  }
  tt <- chron(dates.=sprintf("%04d/%02d/%02d",dts.clm$year,dts.clm$month,dts.clm$day),
              format = c(dates = "y/m/d"))
  dts.clm$doy <- as.numeric(format(as.Date(tt),format="%j"))/365
  dts.clm$anual <- cut(dts.clm$doy,breaks=seq(0,1.01,length=24))

  save(file="~/CEBA/Rdata/JBMts.rda",dts.clm)
}

@ 



%Utilizamos las mediciones diarias de temperatura media, mínima y máxima (°C), humedad relativa media (\%) y precipitación total de lluvia (mm).
 
%Utilizamos datos de climatologías interpoladas para generar una predicción de la temperatura y precipitación esperada en la zona según los datos climáticos de una red de estaciones meteorológicas entre 1960 y 1990 y compararlas con los datos medidos en la estaciones climáticas del área de estudio.

\subsection{Datos de sensores remotos}
Utilizamos productos generados a partir de los sensores Modis Terra y Aqua y disponibles en \citet{MODIS} para caracterizar la temperatura de la superficie de la tierra, índice de vegetación, índice de área foliar y área fotosintéticamente activa, emisividad y albedo, evapotranspiración, etc. Con estos datos generamos series de tiempo para el periodo de 2000 a 2011.

<<DatosModis,echo=false>>=
if (!exists("dts.sen")) {

  dts.sen <- data.frame()
  for (vv in c("250m_16_days_NDVI","250m_16_days_EVI","LST_Day_1km","LST_Night_1km","PET_1km","ET_1km","Lai_1km","Fpar_1km")) {
    for (k in dir(sprintf("%s/%s/%s",mptc,mapaDB,vv))) {
      fch <- sub("A","",strsplit(k,"\\.")[[1]][2])
      yr <- as.numeric(substr( fch,1,4))
      dd <- as.numeric(substr( fch,5,8))
      fch <- yr + (dd/365)
      rq <- raster(sprintf("%s/%s/%s/%s",mptc,mapaDB,vv,k))
      qry <- unlist(extract(rq,JBM))
      dts.sen <- rbind(dts.sen,data.frame(fch=fch,
                                          year=yr,
                                          doy=dd/365,
                                          j=1:length(qry),var=vv,
                                          val=qry))
    }
  }
 
  dts.sen$anual <- cut(dts.sen$doy,breaks=seq(0,1.01,length=24))
  save(file="~/CEBA/Rdata/JBMts.rda",dts.sen,dts.clm)
}
                                                                      
ss <- dts.sen$var %in% c("LST_Day_1km","LST_Night_1km")
dts.sen$val[ss & dts.sen$val < 7500] <- NA
dts.sen$val[ss] <- (dts.sen$val[ss] * 0.02)-273.15
ss <- dts.sen$var %in% c("ET_1km","PET_1km")
dts.sen$val[ss & dts.sen$val > 32760] <- NA
dts.sen$val[ss] <- dts.sen$val[ss] * 0.1
ss <- dts.sen$var %in% c("Fpar_1km","Lai_1km")
dts.sen$val[ss & dts.sen$val > 100] <- NA
dts.sen$val[ss & dts.sen$val < 0] <- NA
dts.sen$val[ss] <- dts.sen$val[ss] * 0.1
ss <- dts.sen$var %in% c("250m_16_days_EVI","250m_16_days_NDVI")
dts.sen$val[ss & dts.sen$val > 10000] <- NA
dts.sen$val[ss & dts.sen$val < -2000] <- NA
dts.sen$val[ss] <- dts.sen$val[ss] * 0.0001
dts.sen$val[ss & dts.sen$val == 0] <- NA

@ 

Enfocamos la comparación en diez pixeles que corresponden con el parche de bosque seco, y diez pixeles representativos de las áreas aledañas con distinto tipo de uso de la tierra.

%Calculamos la cobertura boscosa del fragmento a partir de mapas de un mosaico de imágenes LandSat Hansen et al 2013

\section{Resultados y Discusión}

\subsection{Tendencias climáticas entre 1973 y 2013}

<<SerieTemperatura,echo=false>>=
ss <- subset(dts.clm,year>1972)
ss$fch <- ss$year+ss$doy
@ 

La temperatura media diaria muestra un aparente aumento desde 28° en 1973 hasta 29.7° en 2013, sin embargo las fluctuaciones interanuales son muy marcadas. La temperatura media diaria para todo el periodo es de \Sexpr{sprintf("%0.2f",mean(ss$tmean,na.rm=T))}°C, con valores extremos de \Sexpr{sprintf("%0.2f",max(ss$tmean,na.rm=T))}°C en \Sexpr{with(ss[which.max(ss$tmean),],format.Date(sprintf("%s-%s-%s",year,month,day),"%B de %Y"))} y \Sexpr{sprintf("%0.2f",min(ss$tmean,na.rm=T))}°C en \Sexpr{with(ss[which.min(ss$tmean),],format.Date(sprintf("%s-%s-%s",year,month,day),"%B de %Y"))}.

<<AumentoTemperatura,fig=true,echo=false>>=

##plot(tmean~fch,ss,cex=.25)
##lines(supsmu(ss$fch,ss$tmax,span=.20),col=2,lwd=3)
##lines(supsmu(ss$fch,ss$tmean,span=.20),col=1,lwd=3)
##lines(supsmu(ss$fch,ss$tmin,span=.20),col=4,lwd=3)
##decompose(ts(ss$tmean,start=150,frequency=365))

plot(tmean~fch,ss,cex=.5,pch=1,col="grey55",xlab="Año",ylab="Temperatura [°C]")
ys <- ts(ss$tmean,start=152,frequency=365)
ys <- zoo::na.approx(ys)
lines(ss$fch,decompose(ys)$trend[1:nrow(ss)],col=1,lwd=2)

ys <- ts(ss$tmin,start=152,frequency=365)
ys <- zoo::na.approx(ys)
lines(ss$fch,decompose(ys)$trend[1:nrow(ss)],col=4,lwd=2)

ys <- ts(ss$tmax,start=152,frequency=365)
ys <- zoo::na.approx(ys)
lines(ss$fch,decompose(ys)$trend[1:nrow(ss)],col=2,lwd=2)

@  

Los promedios de las temperaturas máximas y mínimas muestran fluctuaciones similares, sin una tendencia marcada.  El valor promedio de la temperatura máxima diaria es de \Sexpr{sprintf("%0.2f",mean(ss$tmax,na.rm=T))}°C, y la máxima temperatura registrada en todo el período fue de \Sexpr{sprintf("%0.2f",max(ss$tmax,na.rm=T))}°C en \Sexpr{with(ss[which.max(ss$tmax),],format.Date(sprintf("%s-%s-%s",year,month,day),"%B de %Y"))}. En un total de \Sexpr{sum(ss$tmax>40,na.rm=T)} ocasiones la temperatura máxima ha estado por encima de 40°C. El valor promedio de la temperatura mínima diaria es de \Sexpr{sprintf("%0.2f",mean(ss$tmin,na.rm=T))}°C, pero la mínima temperatura registrada fue de \Sexpr{sprintf("%0.2f",min(ss$tmin,na.rm=T))}°C en \Sexpr{with(ss[which.min(ss$tmin),],format.Date(sprintf("%s-%s-%s",year,month,day),"%B de %Y"))}, este último valor fue muy probablemente un error de medición, pues es el único registro en toda la serie por debajo de 10°C, aunque existen al menos  \Sexpr{sum(ss$tmin<15,na.rm=T)} registros de temperatura por debajo de 15°C. % En el caso de la temperatura máxima promedio, el mayor aumento fue entre 1973 y 1989, llegando hasta casi 34°, mientra que la temperatura mínima llegó a un máximo de 25° en el 2000.

<<echo=false>>=
pp <- with(ss,
           aggregate(data.frame(PP=PP),list(yr=year),
                     function(x) {sum(x,na.rm=T)*365*1/length(x)}))
@ 

  La precipitación media anual es de \Sexpr{sprintf("%0.1f",mean(pp$PP))}, pero con una variación interanual de \Sexpr{sprintf("%0.1f",sd(pp$PP))} mm y valores extremos de \Sexpr{sprintf("%0.1f",min(pp$PP))} mm registrado en el año \Sexpr{pp$yr[which.min(pp$PP)]} y \Sexpr{sprintf("%0.1f",max(pp$PP))} mm, en el año \Sexpr{pp$yr[which.max(pp$PP)]}. La media ponderada se encuentra entre 100 y 300 mm en el periodo de 1973 a 1989, y muestra un aumento considerable hasta 400mm anuales entre 1990 y 2009, aunque los valores de los últimos cinco años parecen indicar un descenso a los niveles previos.

<<echo=false,fig=true>>=

plot(PP~yr,pp,type="h",xlab="Año",ylab="Precipitación total anual [mm]")
lines(supsmu(pp$yr,pp$PP),col=3)


@ 

\subsection{Tendencias en vegetación y clima entre 2000 y 2011}

La serie de tiempo del índice de vegetación derivada del satélite MODIS muestra una tendencia positiva en los últimos 11 años a pesar de la marcada estacionalidad. El promedio ponderado se ubica alrededor de 0.5 entre 2000 y 2004 y luego aumenta a valores cercanos o superiores a 0.6. Aparentemente este aumento se debe a un aumento de los valores mínimos anuales con poca variación en los valores máximos, lo cual implica una disminución en la variabilidad intra-anual.

<<NDVIts,echo=false,fig=true>>=
ss <- subset(dts.sen,var %in% "250m_16_days_NDVI")
plot(val~fch,ss,cex=.5,pch=1,xlab="Año",ylab="NDVI")
for (k in 1:46) {
  x <- subset(dts.sen,var %in% "250m_16_days_NDVI" & j ==k)$fch
  y <- subset(dts.sen,var %in% "250m_16_days_NDVI" & j ==k)$val
  ys <- ts(y,start=4,frequency=23)
  ys <- zoo::na.approx(ys)

  lines(x,decompose(ys)$trend,col=2,lwd=2)
}
@ 

Las series de tiempo de la temperatura de la superficie del suelo muestra valores estables de temperatura nocturna entre 23 y 24°C, y una tendencia negativa en la temperatura diurna que disminuye desde 35°C en los años 2000 a 2002 hasta 31° entre 2010 y 2012. 

<<LSTts,echo=false,fig=true>>=
ss <- subset(dts.sen,var %in% c("LST_Day_1km","LST_Night_1km"))
plot(val~fch,ss,cex=.5,pch=1,xlab="Año",ylab="LST")
for (k in 1:46) {
  x <- subset(dts.sen,var %in% "LST_Day_1km" & j ==k)$fch
  if (length(x)>0) {
    y <- subset(dts.sen,var %in% "LST_Day_1km" & j ==k)$val
    ##  x <- x[!is.na(y)]
    ##  y <- y[!is.na(y)]
    ys <- ts(y,start=c(2000,9),frequency=46)
    ys <- zoo::na.approx(ys)
    
    lines(as.numeric(time(ys)),decompose(ys)$trend,col=2,lwd=2)
  }
  x <- subset(dts.sen,var %in% "LST_Night_1km" & j ==k)$fch
  if (length(x)>0) {
    
    y <- subset(dts.sen,var %in% "LST_Night_1km" & j ==k)$val
    ##  x <- x[!is.na(y)]
    ##  y <- y[!is.na(y)]
    ys <- ts(y,start=c(2000,9),frequency=46)
    ys <- zoo::na.approx(ys)
    
    lines(as.numeric(time(ys)),decompose(ys)$trend,col=4,lwd=2)
  }
}
@ 

\subsection{Estacionalidad en el periodo 2000 a 2011}


Ver \citet{Vila60} para los detalles del patrón anual de precipitación.

%El patrón de precipitación anual es similar al periodo 

Entre 2000 y 2011 los promedios quincenales de temperatura máxima y mínima en la estación climática muestran poca estacionalidad y poca variación interanual con valores ligeramente más altos entre julio y septiembre.

<<temperatura,echo=false,fig=true>>=
##tt <- with(subset(dts.clm,year>1972),aggregate(data.frame(tmean=tmean,tmin=tmin,tmax=tmax),list(yr=year,anual=anual),function(x) {mean(x,na.rm=T)}))

tt <- with(subset(dts.clm,year>1999 & year<2012),
           aggregate(data.frame(tmean=tmean,tmin=tmin,tmax=tmax),
                     list(yr=year,anual=anual),function(x) {mean(x,na.rm=T)}))
tx <- with(subset(dts.clm,year>1999 & year<2012),
           aggregate(data.frame(tmax=tmax),list(yr=year,anual=anual),
                     function(x) {max(x,na.rm=T)}))
tn <- with(subset(dts.clm,year>1999 & year<2012),
           aggregate(data.frame(tmin=tmin),list(yr=year,anual=anual),
                     function(x) {min(x,na.rm=T)}))

par(xpd=T)
boxplot(tmax~anual,tx,ylim=c(10,45),axes=F,ylab="Temperatura [°C]")
boxplot(tmin~anual,tn,axes=F,add=T)
boxplot(tmean~anual,tt,ylim=c(20,35),axes=F,add=T)
axis(2)
xxs <- 0.5+(23*c(0,31,59,90,120,151,181,212,243,273,304,334,365)/365)
for (j in 1:12)
  polygon(c(xxs[j],xxs[j],xxs[j+1],xxs[j+1],xxs[j]),
          c(-1.65,-.5,-.5,-1.65,-1.65)+10)
text(xxs[1:12]+(23*16/365),9,c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dec"))

@ 


<<echo=false,results=hide>>=
summary(with(subset(dts.clm,year>1972),aggregate(data.frame(tmean=tmean,tmin=tmin,tmax=tmax),list(yr=year),function(x) {max(x,na.rm=T)})))
summary(with(subset(dts.clm,year>1972),aggregate(data.frame(tmean=tmean,tmin=tmin,tmax=tmax),list(yr=year),function(x) {min(x,na.rm=T)})))
@ 



%Las mediciones de temperatura de la estación climática muestran un aumento estacional muy regular en los meses de … a... para la temperatura media, máxima y mínima, pero los estimados de temperatura en la superficie de la tierra basadas en sensores remotos muestran gran variabilidad intranual en la temperatura diurna y más regularidad en la nocturna. De hecho, a pesar de que el valor de LST nocturno se corresponden muy bien con las mediciones de temperatura mínima, el valor de LST diurno muestra picos en los meses de marzo seguido por valores bajos entre … y ….


%\subsection{Precipitacion E.M.}


<<echo=false,results=hide>>=
ss <- subset(dts.clm,year>1999 & year<2012) ##subset(dts.clm,year>1972),
pttl <- with(ss,aggregate(data.frame(PP=PP),list(yr=year),function(x) {sum(x,na.rm=T)*365/length(x)}))
pver <- with(subset(ss, anual %in% levels(anual)[7:13]),
             aggregate(data.frame(PP=PP),list(yr=year),function(x) {sum(x,na.rm=T)*365/23*7/length(x)}))
phum <- with(subset(ss, anual %in% levels(anual)[14:21]),
             aggregate(data.frame(PP=PP),list(yr=year),function(x) {sum(x,na.rm=T)*365/23*8/length(x)}))

pttl <- merge(pttl,merge(phum,pver,by="yr"),by="yr")
summary(pttl$PP.x)
summary(pttl$PP.x/pttl$PP)
summary(pttl$PP.y/pttl$PP)
summary((pttl$PP.x+pttl$PP.y)/pttl$PP)

@ 
La precipitación quincenal si muestra un patrón estacional marcado, con un periodo de escasas precipitaciones entre la segunda mitad de diciembre y el comienzo de abril. Se observa un primer periodo de lluvias entre mayo y junio, y luego un pico más pronunciado entre septiembre y octubre. En promedio, el \Sexpr{sprintf("%0.1f",mean(pttl$PP.x/pttl$PP)*100)}\% de la precipitación anual cae entre agosto y noviembre, y el \Sexpr{sprintf("%0.1f",mean(pttl$PP.y/pttl$PP)*100)}\% cae entre abril y julio, sin embargo estos porcentajes pueden variar, y en algunos años (notablemente en \Sexpr{pttl$yr[with(pttl,which.max(PP.y/PP))]}), el periodo de abril a julio aporta mayor precipitación que el resto del año. 
  
<<PrecipitacionTotalAnual,echo=false,fig=true>>=

pp <- with(ss,
           aggregate(data.frame(PP=PP),list(yr=year,anual=anual),
                     function(x) {sum(x,na.rm=T)*365/23*1/length(x)}))

par(xpd=T)
boxplot(PP~anual,pp,ylim=c(0,200),axes=F,ylab="Precipitación [mm]")
axis(2)
xxs <- 0.5+(23*c(0,31,59,90,120,151,181,212,243,273,304,334,365)/365)
for (j in 1:12)
  polygon(c(xxs[j],xxs[j],xxs[j+1],xxs[j+1],xxs[j]),
          c(-2.15,-.5,-.5,-2.15,-2.15)*5)
text(xxs[1:12]+(23*16/365),-7,c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dec"))
box(,bty="n")

@ 




<<echo=false,eval=false>>=
##dias con precipitacion mayor a 1mm
aggregate(dts.clm$PP>1,list(dts.clm$year),sum,na.rm=T)

@ 
%\subsection{Temperatura Modis}

En contraste con la temperatura medida en la estación climática, los estimados de temperatura en la superficie terrrestre (LST) medidos con sensores remotos muestran una estacionalidad mucho más marcada. Las mayores temperaturas diurnas se observan entre marzo y mayo, y fluctua en los próximos meses hasta llegar a valores mínimos entre octubre y noviembre. La variación interanual es muy marcada, especialmente en abril y entre agosto y septiembre. Por su parte la temperatura nocturna muestra un patrón unimodal con un máximo entre agosto y septiembre, similar a las mediciones realizadas en la estación meteorológica.

<<Tmodis,echo=false,fig=true>>=
tx <- with(subset(dts.sen,var %in% "LST_Day_1km"),aggregate(data.frame(tmax=val),list(yr=year,anual=anual),function(x) {max(x,na.rm=T)}))
tn <- with(subset(dts.sen,var %in% "LST_Night_1km"),aggregate(data.frame(tmin=val),list(yr=year,anual=anual),function(x) {min(x,na.rm=T)}))

tx <- subset(tx,is.finite(tx$tmax))
tn <- subset(tn,is.finite(tn$tmin))
par(xpd=T)
boxplot(tmax~anual,tx,ylim=c(10,45),axes=F,ylab="Temperatura [°C]",
        col="grey77",border="grey35")
boxplot(tmin~anual,tn,axes=F,add=T,
        col="grey12",border="grey35")
##boxplot(tmax~anual,tx,ylim=c(20,35),axes=F,add=T)
axis(2)
xxs <- 0.5+(23*c(0,31,59,90,120,151,181,212,243,273,304,334,365)/365)
for (j in 1:12)
  polygon(c(xxs[j],xxs[j],xxs[j+1],xxs[j+1],xxs[j]),
          c(-1.65,-.5,-.5,-1.65,-1.65)+10)
text(xxs[1:12]+(23*16/365),9,c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dec"))

@ 


<<echo=false,results=hide>>=
summary(aggregate(tt$tmean,list(tt$yr),mean,na.rm=T)$x)
summary(aggregate(tx$tmax,list(tx$yr),max,na.rm=T)$x)
summary(aggregate(tn$tmin,list(tn$yr),min,na.rm=T)$x)

@ 
%\subsection{Evapotranspiracion Modis}

Los valores de evapotranspiración estimados a partir de imágenes de satélite muestran dos periodos bien diferenciados. Durante los meses más secos (febrero y marzo) los valores llegan a un valor mínimo de 10, pero luego aumenta y entre junio y noviembre se situa alrededor de 40, aunque con notable variación interanual.

<<ET,echo=false,fig=true>>=
ee <- with(subset(dts.sen,var %in% "ET_1km"),aggregate(data.frame(ET=val),list(yr=year,anual=anual),function(x) {max(x,na.rm=T)}))
ep <- with(subset(dts.sen,var %in% "PET_1km"),aggregate(data.frame(PET=val),list(yr=year,anual=anual),function(x) {max(x,na.rm=T)}))


par(xpd=T)
boxplot(ET~anual,ee,ylim=c(0,100),axes=F,ylab="Evapotranspiración []",
        col="grey77",border="grey35")
##boxplot(PET~anual,ep,axes=F,add=T,
##        col="grey12",border="grey35")
##boxplot(tmax~anual,tx,ylim=c(20,35),axes=F,add=T)
axis(2)
xxs <- 0.5+(23*c(0,31,59,90,120,151,181,212,243,273,304,334,365)/365)
for (j in 1:12)
  polygon(c(xxs[j],xxs[j],xxs[j+1],xxs[j+1],xxs[j]),
          c(-3.65,-.5,-.5,-3.65,-3.65)+0)
text(xxs[1:12]+(23*16/365),-2,c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dec"))

@ 

%\subsection{NDVI Modis}

Por su parte el índice de vegetación sugiere que la vegetación tiene una fenología bimodal ligeramente desfasada de la precipitación, con un mínimo al final de la sequía y comienzo de las primeras lluvias (entre marzo y abril), un primer pico en junio, con valores sostenidos hasta septiembre y un segundo incremento al final de la época de lluvias entre octubre y noviembre.

<<NDVI,echo=false,fig=true>>=
nn <- with(subset(dts.sen,var %in% "250m_16_days_NDVI"),aggregate(data.frame(NDVI=val),list(yr=year,anual=anual),function(x) {mean(x,na.rm=T)}))

par(xpd=T)
boxplot(NDVI~anual,nn,ylim=c(0,1),axes=F,ylab="NDVI",
        col="grey77",border="grey35")
##boxplot(PET~anual,ep,axes=F,add=T,
##        col="grey12",border="grey35")
##boxplot(tmax~anual,tx,ylim=c(20,35),axes=F,add=T)
axis(2)
xxs <- 0.5+(23*c(0,31,59,90,120,151,181,212,243,273,304,334,365)/365)
for (j in 1:12)
  polygon(c(xxs[j],xxs[j],xxs[j+1],xxs[j+1],xxs[j]),
          c(-.1,-.05,-.05,-.1,-.1)+0)
text(xxs[1:12]+(23*16/365),-.075,c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dec"))

@ 

%\section{Correlaciones}

<<echo=false,eval=false>>=
 with(merge(ep,pp,by=c("yr","anual")),cor.test(PP,PET))
with(merge(ee,pp,by=c("yr","anual")),cor.test(PP,ET))

@ 
%La evapo-transpiración sigue un patrón similar a las precipitaciones, con valores mínimos de … y máximos entre … y …
%Sin embargo la evapotranspiración potencial es mucho mayor a principios de año, llegando a estar muy por encima de la evapotranspiración estimada, en el segundo semestre del año ET y ETP tienen valores similares.

%El índice de vegetación tiene valores mínimos en marzo (0.2) y valores máximos entre octubre y noviembre (0.6), Entre mayo y enero el valor de NDVI suele mantenerse por encima de 0.4, aunque se observa una leve reducción a mediados de julio, que coincide con un breve periodo de reducción de la precipitación. 

%El climadiagrama normalizado para el periodo 1960-1990 muestra una temperatura media anual de 27.7°C y una precipitación total anual de 580 mm (215-1250mm), con 4.2 días con precipitación mayor o igual a 1.0 mm, evaporación total de 1826 mm, humedad relativa de 76\% (50-92\%) \citep{NormClim1993}. El valor total de precipitación obtenido a partir de la página web no corresponde exactamente con los totales disponibles en documentos de ...

%Los datos medidos por sensores remotos estiman una temperatura promedio diurna de … y nocturna de …, el valor diurno generalmente está por encima de los valores de temperatura máxima medidos en la estación climática, pero los valores diurnos corresponden muy bien con los valores de temperaturas mínimas. 

%la baja precipitación anual 
%La vegetación usualmente tiene una respuesta moderada a los cambios en las precipitaciones y 
%A pesar de las variaciones intra-anuales en la precipitación,  sin embargo las variaciones en la precipitación entre los meses de mayo a noviembre no son 

%Se observa claramente que la vegetación sigue un perfil uni/bi-modal con picos máximos en … y mínimos en … 
%La tendencia en el índice de vegetación es creciente




%Decreto y reglamento de uso de la Zona Protectora del Área Metropolitana de la ciudad de Maracaibo entre 1986 y 1989\citep{DecretoZPM,ReglamentoZPM}.

%Zona protectora de Maracaibo

%La Zona Protectora de Maracaibo es una franja territorial de unas 20.000 hectáreas que rodea a la ciudad de Maracaibo en Venezuela, a modo de cinturón verde (green-belt) para coadyuvar en la contención de la expansión urbana y la protección de recursos naturales que rodean a la ciudad. Fue decretada por la Presidencia de la República de Venezuela a finales de 1980, y permanece como territorio bajo jurisdicción del Ministerio de Poder Popular para el Ambiente y los Recursos Naturales de Venezuela.

%La Zona Protectora de Maracaibo califica dentro de las denominadas Áreas Bajo Régimen de Administración Especial, en el marco jurídico de la legislación territorial venezolana.


%Esta zona pertenece a la Zona Protectora de Maracaibo,
%En el caso de la Zona Protectora de la ciudad de Maracaibo, este concepto de “cinturón verde” prevalece cuando se plantea que la creación de ZP “tiene como objetivo primordial frenar la expansión anárquica de la ciudad y actuar como agente regulador del clima y el medio ambiente en beneficio de la calidad de la vida urbana. Además contempla la protección y el uso adecuado de los recursos naturales, así como la creación de alternativas de recreación para los habitantes de la ciudad de Maracaibo” (Plan de Ordenación Territorial Zona Protectora Ciudad Maracaibo, 1989: 3). Asimismo, se permiten y propician los usos agrícolas u otros usos complementarios y compatibles con el ambiente natural de la zona.
\section{Discusión}

Se observó una tendencia positiva en la temperatura entre 1970 y el presente, pero en los últimos 10 años una tendencia positiva en la cobertura vegetal con disminución de la temperatura diurna de la superficie de la tierra. A pesar del régimen de precipitaciones bimodal, la vegetación en el área de estudio muestra un pefil intermedio entre un régimen unimodal y uno bimodal con un periodo prolongado de crecimiento vegetativo. El regimen de precipitación es variable entre años, y la vegetación es capaz de mantenerse activa durante varios meses con variaciones de precipitación, pero las condiciones no son suficientes para mantener una cobertura boscosa más densa.  


%%\section*{Referencias}
%% en vez de printbibliography
%%\bibliography{/home/jferrer/CEBA/lib/BibTEX/Suelos}
\bibliography{/home/jferrer/CEBA/lib/BibTEX/JBM}


\end{document}