#Citar: Nuñez, J. 2011. RSARFLM v.1. Regional Frequency Analysis L-moments R Script. Water Center for Arid and Semiarid Zones of Latina America and the Caribbean.CAZALAC.La Serena,Chile #RS-RFA-LM-CAZALAC: VERSION 2.12-14-11-2011 # ANALISIS REGIONAL DE FRECUENCIA DE LA PRECIPITACION BASADO EN L-MOMENTOS #AUTOR: JORGE NUÑEZ COBO, CAZALAC, PROGRAMA DOCTORADO ING.AGR.UNIVERSIDAD DE CONCEPCION-CHILE # REFERENCIAS: Nuñez,J., Verbist,K., Wallis,J., Schaefer,M., Morales,L. and Cornelis,W. 2011.Regional frequency analysis for mapping drought events in north-central Chile.Journal of hydrology. 405 (3-4), pg. 352-366 # REFERENCIAS: UNESCO. 2010. Guía Metodológica para la Aplicación del Análisis Regional de Frecuencia de Sequías basado en L-Momentos y Resultados de Aplicación en América Latina. Documento Técnico PHI-LAC Nº 27. ------------------------------------------------------------------------------------------------------- # BLOQUE 1: PREPARACION DEL SISTEMA #------------------------------------------------------------------------------------------------------- # PASO 1: INSTALAR PAQUETES # Si no se tiene los paquetes instalados, eliminar símbolo de comentario (#) para activar los comandos abajo #install.packages("lmom")# Original de Hosking #install.packages("lmomRFA") # Segundo paquete de Hosking #install.packages("nsRFA") # Alberto Viglione #install.packages("raster") # Trabajar con imagenes raster, para etapa de mapeo #install.packages("rgdal")# Formatos de archivos #install.packages("sp") # Paquete para objetos espaciales #install.packages("DEoptim") # Paquete para minimizacion ajuste curvas #install.packages("sqldf")# Paquete para manipulación de Base de Datos #install.packages("tcltk")# Paquete para seleccionar carpeta de trabajo #install.packages("Matrix")# Paquete para calcular el número de valores no cero #install.packages("maptools")# Paquete para calcular el número de valores no cero #install.packages("reshape")# Paquete para calcular el número de valores no cero #install.packages("ggplot2")# Paquete para calcular el número de valores no cero #install.packages("car") # PASO 2: CARGAR PAQUETES library(lmom) library(lmomRFA) library(nsRFA) library(raster) library(rgdal) library(sp) library(DEoptim) library(sqldf) library(tcltk) library(Matrix) library(maptools) library(reshape) library(ggplot2) library(car) #NOTA: SI LA BASE DE DATOS YA EXISTE Y ESTA DENTRO DE LA CARPETA DE TRABAJO, EJECUTAR LAS LINEAS 39-40 Y SEGUIR EN LA LINEA 143 # Los requisitos para ejecutar todo el procedimiento son: # a) Un archivo de base de datos (cvs) con registros de la variable de interés # b) Un archivo de base de datos (csv) con características de las estaciones a utilizar # c) Un mapa en formato raster (geotiff) de la variable explicativa (precipitación media anual) #PASO 3: SELECCION DE LA CARPETA DE TRABAJO DONDE SE ALMACENARA TODA LA INFORMACIÓN (archivos,mapas) WF<-tk_choose.dir(getwd(), "Choose a suitable folder")#Carpeta de trabajo lo más cerca de MIS DOCUMENTOS setwd(WF) ------------------------------------------------------------------------------------------------------- # BLOQUE 2: ESTRUCTURACION DE LA BASE DE DATOS #------------------------------------------------------------------------------------------------------- # EL PRESENTE SCRIPT UTILIZA UNA BASE DE DATOS ESTRUCTURADA EN DOS ARCHIVOS: UNO CON REGISTROS Y OTRO CON # CARACTERISTICAS DE LAS ESTACIONES # PASO 4: IMPORTACION DE LA BASE DE DATOS INICIAL AL SISTEMA Y ASEGURAMIENTO DE CALIDAD # Primero leo la base de datos que está constituida por dos archivos separados. # La primera almacena información de los registros y la segunda, sobre las estaciones. ISO_PAIS<-"ALC" # En este comando, entre paréntesis, se indica el código ISO de 3 dígitos del país según normal ISO 3166-1 (http://es.wikipedia.org/wiki/ISO_3166-1) # El código ISO sera incluido en la salida de los archivos. Digitar en MAYUSCULAS. BaseDatosEstaciones <- read.csv("BaseDatosEstaciones.csv", sep=";",na.strings="NA") #IMPORTANTE: EL NOMBRE DEL ARCHIVO DE ENTRADA DEBE SER EXACTAMENTE BaseDatosEstaciones.csv BaseDatosRegistros <- read.csv("BaseDatosRegistros.csv", sep=";",na.strings="NA")#IMPORTANTE: EL NOMBRE DEL ARCHIVO DE ENTRADA DEBE SER EXACTAMENTE BaseDatosRegistros.csv #ANALIZAR SI LAS BASES DE DATOS ESTAN OK Y SIN ERROR EN SU ESTRUCTURA str(BaseDatosEstaciones) # Analizar que el tipo de valor que hay en cada columna sea el que corresponde. # Si debe haber número en una columna pero indica que hay "factor" es que seguramente se incluyó texto en vez de número en la columna correspondiente str(BaseDatosRegistros)# Realiza el mismo análisis, pero para la base de Registros # SI NO, CONVERTIR CUALQUIER DATO QUE NO ES NUMERO, EN NA BaseDatosRegistros$anio<- as.numeric(as.character(BaseDatosRegistros$anio)) BaseDatosRegistros$ENE <- as.numeric(as.character(BaseDatosRegistros$ENE)) BaseDatosRegistros$FEB <- as.numeric(as.character(BaseDatosRegistros$FEB)) BaseDatosRegistros$MAR <- as.numeric(as.character(BaseDatosRegistros$MAR)) BaseDatosRegistros$ABR <- as.numeric(as.character(BaseDatosRegistros$ABR)) BaseDatosRegistros$MAY <- as.numeric(as.character(BaseDatosRegistros$MAY)) BaseDatosRegistros$JUN <- as.numeric(as.character(BaseDatosRegistros$JUN)) BaseDatosRegistros$JUL <- as.numeric(as.character(BaseDatosRegistros$JUL)) BaseDatosRegistros$AGO <- as.numeric(as.character(BaseDatosRegistros$AGO)) BaseDatosRegistros$SEP <- as.numeric(as.character(BaseDatosRegistros$SEP)) BaseDatosRegistros$OCT <- as.numeric(as.character(BaseDatosRegistros$OCT)) BaseDatosRegistros$NOV <- as.numeric(as.character(BaseDatosRegistros$NOV)) BaseDatosRegistros$DIC <- as.numeric(as.character(BaseDatosRegistros$DIC)) # PASO 5: CALCULO DE VARIABLES DERIVADAS Y CONSTRUCCION DE LA BASE DE DATOS DEFINITIVA # A continuación genero la base de datos útil, constituida por los registros que no contienen datos nulos. # Además, determino la cantidad de estaciones útiles EstacionesOriginales<-as.factor(BaseDatosRegistros[[1]]) NumeroEstacionesOriginales<-nlevels(EstacionesOriginales)# Cada estación es 1 nivel PP<-na.omit(BaseDatosRegistros, row.names=F)# Me quedo solamente con registros completos. EstacionesCompletas<-as.factor(PP[[1]]) NumeroEstacionesCompletas<-nlevels(EstacionesCompletas)# Comparo este valor con "6" para ver efecto sobre tamaño de la Base de Datos # Acá leo las estaciones con sus códigos originales. LluviaAnual<-PP[3:14] # Calculo la lluvia anual. Pero podria usar otra variable y cambiar las columnas. # Ejemplo: LluviaJulio<-PPALC[9] # Ejemplo: LluviaInvierno<-PPALC[7:10] L<-length(PP[[1]]) # Obtengo el valor de longitud total de registros SumaLluviaAnual<-matrix(rowSums(LluviaAnual[1:L,]),nrow=L,ncol=12) SumaLluviaAnual<-ifelse(SumaLluviaAnual>0,SumaLluviaAnual,0.1)# # ..........INICIO CALCULO DE INDICE DE ESTACIONALIDAD (SI) Y DIA JULIANO MEDIO (JMD).............. # Aca se inicia el calculo del Indice de Estacionalidad y Dia Juliano. DiaJulianoAng<-matrix(seq(15,345,30)*2*pi/365,nrow=L,ncol=12,byrow=TRUE) Prec<-PP[1:L,3:14] x<-Prec*cos(DiaJulianoAng)/SumaLluviaAnual y<-Prec*sin(DiaJulianoAng)/SumaLluviaAnual xcos<-matrix(rowMeans(x),nrow=L,ncol=1) ysin<-matrix(rowMeans(y),nrow=L,ncol=1) xcossum<-matrix(rowSums(x),nrow=L,ncol=1) ysinsum<-matrix(rowSums(y),nrow=L,ncol=1) angulo<-atan(ysin/xcos) angulo_corregido<-matrix(0,nrow=L,ncol=1) for (k in 1:L) { if (xcos[k]>0&ysin[k]>0) angulo_corregido[k]<-angulo[k] else if (xcos[k]<0&ysin[k]<0) angulo_corregido[k]<-angulo[k]+pi else if (ysin[k]>0&xcos[k]<0) angulo_corregido[k]<-angulo[k]+pi else angulo_corregido[k]<-angulo[k]+2*pi } JMD<-(angulo_corregido*365)/(2*pi) SI<-sqrt(xcossum^2+ysinsum^2) # ............ FIN DE CALCULO INDICE DE ESTACIONALIDAD Y DIA JULIANO MEDIO........................ BaseDatosIntermedia<-cbind(PP,SumaLluviaAnual[,1],SI,JMD)# Uno las columnas con los Indices calculados para cada registro colnames(BaseDatosIntermedia)[15]<-'SumaLluviaAnual' # Acá inicio el cálculo de los valores medios de cada estación #Tengo que usar tapply para calcular las medias de SI y MAP y JMD y luego estos valores asignarlos a las estaciones #y esa matriz la uso para regionalizar según SI, MAP o JMD. # Esta calcula el SI medio por estación.Sin datos ausentes. SI_por_Estacion<-as.matrix(tapply(BaseDatosIntermedia[[16]],BaseDatosIntermedia[[1]],mean,na.rm=TRUE)) hist(SI_por_Estacion)# Visualización de resultados para detectar fallas de cálculo de SI # Esta calcula la PMA medio por estación.Sin datos ausentes. PMA_por_Estacion<-as.matrix(tapply(BaseDatosIntermedia[[15]],BaseDatosIntermedia[[1]],mean,na.rm=TRUE)) hist(PMA_por_Estacion)# Visualización de resultados para detectar fallas de calculo de MAP # Esta calcula el JMD medio por estación.Sin datos ausentes. JMD_por_Estacion<-as.matrix(tapply(BaseDatosIntermedia[[17]],BaseDatosIntermedia[[1]],mean,na.rm=TRUE)) hist(JMD_por_Estacion)# Visualización de resultados para detectar fallas de cálculo de JMD # Esta calcula la Longitud de Registro media por estación.Sin datos ausentes. LR_por_Estacion<-as.matrix(tapply(BaseDatosIntermedia[[15]],BaseDatosIntermedia[[1]],length)) hist(LR_por_Estacion)# Visualización de resultados para detectar fallas de cálculo de RL # Esta calcula la cantidad de datos no cero media por estación.Sin datos ausentes. NonZeroNumbers_por_Estacion_Enero<-as.matrix(tapply(BaseDatosIntermedia[[3]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosENE<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Enero)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Febrero<-as.matrix(tapply(BaseDatosIntermedia[[4]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosFEB<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Febrero)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Marzo<-as.matrix(tapply(BaseDatosIntermedia[[5]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosMAR<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Marzo)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Abril<-as.matrix(tapply(BaseDatosIntermedia[[6]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosABR<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Abril)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Mayo<-as.matrix(tapply(BaseDatosIntermedia[[7]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosMAY<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Mayo)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Junio<-as.matrix(tapply(BaseDatosIntermedia[[8]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosJUN<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Junio)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Julio<-as.matrix(tapply(BaseDatosIntermedia[[9]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosJUL<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Julio)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Agosto<-as.matrix(tapply(BaseDatosIntermedia[[10]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosAGO<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Agosto)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Septiembre<-as.matrix(tapply(BaseDatosIntermedia[[11]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosSEP<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Septiembre)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Octubre<-as.matrix(tapply(BaseDatosIntermedia[[12]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosOCT<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Octubre)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Noviembre<-as.matrix(tapply(BaseDatosIntermedia[[13]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosNOV<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Noviembre)/LR_por_Estacion # IMPORTANTE NonZeroNumbers_por_Estacion_Diciembre<-as.matrix(tapply(BaseDatosIntermedia[[14]],BaseDatosIntermedia[[1]],nnzero)) #IMPORTANTE PropCerosDIC<-(LR_por_Estacion-NonZeroNumbers_por_Estacion_Diciembre)/LR_por_Estacion # IMPORTANTE BaseDatosCeros<-cbind(PropCerosENE,PropCerosFEB,PropCerosMAR,PropCerosABR,PropCerosMAY,PropCerosJUN,PropCerosJUL,PropCerosAGO,PropCerosSEP,PropCerosOCT,PropCerosNOV,PropCerosDIC) colnames(BaseDatosCeros)[c(1,2,3,4,5,6,7,8,9,10,11,12)]<-c('PROPCEROSENE','PROPCEROSFEB','PROPCEROSMAR','PROPCEROSABR','PROPCEROSMAY','PROPCEROSJUN','PROPCEROSJUL','PROPCEROSAGO','PROPCEROSSEP','PROPCEROSOCT','PROPCEROSNOV','PROPCEROSDIC') #write.csv(BaseDatosCeros, file = "BaseDatosCeros.csv",row.names=TRUE) id_estacion<-levels(EstacionesCompletas)# Identifico las estaciones a utilizar BaseDatosIndices<-cbind(id_estacion,SI_por_Estacion,PMA_por_Estacion,JMD_por_Estacion,LR_por_Estacion, BaseDatosCeros)# Unos los valores calculados para cada estación colnames(BaseDatosIndices)[c(2,3,4,5)]<-c('SIMedio','PMA','JMDMedio','LongRegistro') # Asigno los nombres de las nuevas columnas calculadas BaseConsolidada<-merge(BaseDatosEstaciones,BaseDatosIndices,by.x="id_estacion",by.y="id_estacion") # Acá uno las dos bases de datos a nivel de estaciones solamente. BaseConsolidada_sin_NA<-na.omit(BaseConsolidada) # Elimino todas las estaciones con datos ausentes. write.csv(BaseConsolidada_sin_NA, file = paste("BaseConsolidada",ISO_PAIS,"sin_NA.csv",sep=''),row.names=FALSE) #Para crear una sola gran Base de Datos BaseCompleta<-merge(BaseConsolidada_sin_NA,BaseDatosIntermedia, by.x = "id_estacion", by.y = "id_estacion") # Acá uno la Base de Datos de Estaciones con Indices Medios y la Base de Datos de Registros write.csv(BaseCompleta, file = paste("BaseCompleta",ISO_PAIS,".csv",sep=''),row.names=FALSE) #---------------------------------------------------------------------------------------------------------------------------- # BLOQUE 3: CREACION DE LAS REGIONES HOMOGENEAS "A PRIORI" #---------------------------------------------------------------------------------------------------------------------------- # PASO 6: REACTUALIZACION DE LA BASE DE DATOS remove(BaseCompleta) # Aca debo hacer esto para reactualizar la Base de Datos en el sistema remove(BaseConsolidada_sin_NA) # Si la Base de Datos ya existe, partir acá BaseCompleta <- read.csv(file=paste("BaseCompleta",ISO_PAIS,".csv",sep=''))# Y vuelvo a cargar la Base de Datos BaseConsolidada_sin_NA<-read.csv(paste("BaseConsolidada",ISO_PAIS,"sin_NA.csv",sep='')) # PASO 7: CREACION DE LAS REGIONES HOMOGENEAS SEGUN DISTINTOS CRITERIOS MEDIANTE COMANDOS SQL CON PAQUETE SQLDF #Acá empiezo a crear las regiones según distintos criterios # La recomendación será: Agrupar por SIMedio (0-0.2, 0.2-0.4, 0.4-0.6,0.6-0.8,0.8-1) # Luego, en cada grupo de SIMedio, separar por JMD (grupos de 30 dias) # Y dentro de este grupo, separar por rangos de Precipitación Media Anual (MAP) # CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC # Aca preparo una región TOTAL para analizar toda ALC de una vez. RegionTotal<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompleta where LongRegistro>14 ") RegionTotal_dat<-RegionTotal["SumaLluviaAnual"][,] RegionTotal_fac<-factor(RegionTotal["id_estacion"][,]) RegTotal<-split(RegionTotal_dat,RegionTotal_fac)# Con esto separo los registros según la estación # CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC # Para autocorrelacion #durbin.watson {car} #............................................................................................ # Acá hago las regiones una por una Region1<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 50 and 159 and LR>15") Region1_dat<-Region1["SumaLluviaAnual"][,] Region1_fac<-factor(Region1["id_estacion"][,]) Reg1<-split(Region1_dat,Region1_fac)# Con esto separo los registros según la estación Region2<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 160 and 227 and LR>15") Region2_dat<-Region2["SumaLluviaAnual"][,] Region2_fac<-factor(Region2["id_estacion"][,]) Reg2<-split(Region2_dat,Region2_fac) Region3<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 227 and 261 and LR>15") Region3_dat<-Region3["SumaLluviaAnual"][,] Region3_fac<-factor(Region3["id_estacion"][,]) Reg3<-split(Region3_dat,Region3_fac) Region4<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 261 and 306 and LR>15") Region4_dat<-Region4["SumaLluviaAnual"][,] Region4_fac<-factor(Region4["id_estacion"][,]) Reg4<-split(Region4_dat,Region4_fac) Region5<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 306 and 396 and LR>15") Region5_dat<-Region5["SumaLluviaAnual"][,] Region5_fac<-factor(Region5["id_estacion"][,]) Reg5<-split(Region5_dat,Region5_fac) Region6<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 396 and 463 and LR>15") Region6_dat<-Region6["SumaLluviaAnual"][,] Region6_fac<-factor(Region6["id_estacion"][,]) Reg6<-split(Region6_dat,Region6_fac) Region7<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 463 and 566 and LR>15") Region7_dat<-Region7["SumaLluviaAnual"][,] Region7_fac<-factor(Region7["id_estacion"][,]) Reg7<-split(Region7_dat,Region7_fac) Region8<-sqldf("select id_estacion, SumaLluviaAnual from BaseCompletaALC where PMA between 566 and 1215 and LR>15") Region8_dat<-Region8["SumaLluviaAnual"][,] Region8_fac<-factor(Region8["id_estacion"][,]) Reg8<-split(Region8_dat,Region8_fac) # ___________________________________________________________________________________________ # Generar algunas graficas interesantes de distribución mensual para una estación de interés #OPCION 1: EN ESTA OPCION GENERO GRAFICAS PARA UN CONJUNTO DE ESTACIONES SELECCIONADAS SEGUN CRITERIOS QUE UNO DEFINE. # ES DECIR, SELECCIONO LAS ESTACIONES QUE ME INTERESAN,PARA ANALIZAR GRAFICAMENTE UN PATRON ListadoEstaciones<-sqldf("select nombre_estacion from BaseCompleta where LongRegistro>50")# ACA SELECCIONAR LOS CRITERIOS PARA SELECCIONAR LAS ESTACIONES LEst<-levels(factor(ListadoEstaciones["nombre_estacion"][,])) Listado<-as.character(LEst) par(mfrow=c(1,2)) for (m in 1:length(Listado)){ Boxplotmensual<-sqldf(paste("select ENE,FEB,MAR,ABR,MAY,JUN,JUL,AGO,SEP,OCT,NOV,DIC from BaseCompleta where nombre_estacion=='",Listado[m],"'",sep="")) boxplot(Boxplotmensual, main=Listado[m],pch=20,cex=0.8, cex.main=1, cex.axis=0.5, cex.lab=0.9, las=1, xlab="Mes", ylab="Precipitación[mm]") BP<-matrix(colMeans(Boxplotmensual),nrow=1) barplot(BP,names.arg=c("E","F","M","A","M","J","J","A","S","O","N","D"),main=Listado[m], cex.names=0.5,cex.main=1, xlab="Mes", ylab="Precipitación[mm]") } plot.new() #OPCION 2: SELECCIONO UNA ESTACION ESPECIFICA A GUSTO PROPIO. SI SE REQUIEREN MAS DE UNA, COPIAR Y PAGAR ABAJO Y CAMBIAR NOMBRE DE ESTACION par(mfrow=c(1,2)) Boxplotmensual<-sqldf("select ENE,FEB,MAR,ABR,MAY,JUN,JUL,AGO,SEP,OCT,NOV,DIC from BaseCompleta where nombre_estacion=='OVALLE'") boxplot(Boxplotmensual, main="St. Ovalle",pch=20,cex=0.8, cex.main=1, cex.axis=0.5, cex.lab=0.9, las=1, xlab="Mes", ylab="Precipitación[mm]") BP<-matrix(colMeans(Boxplotmensual),nrow=1) barplot(BP,names.arg=c("E","F","M","A","M","J","J","A","S","O","N","D"),main="St. Ovalle", cex.names=0.5,cex.main=1, xlab="Mes", ylab="Precipitación[mm]") plot.new() # Ejemplo para seleccionar una estación en particular #RegionXX <- sqldf("select * from BaseCompleta where id_estacion=='st-nnn-0001'") # Ejemplo para seleccionar todas menos una determinada estación #Regionzz<- sqldf("select * from BaseCompleta where id_estacion!='st-nnn-0001'") # Referencia: Halekoh et al, 2010. Handling large(r) datasets in R. http://genetics.agrsci.dk/~sorenh/misc/Rdocs/R-largedata.pdf # otra opción es usando las regiones ingresadas entre líneas 30-46 #_____________________________________________________________________________________________________________ BaseRegiones<-list(RegTotal)#Reg1,Reg2,Reg3,Reg4,Reg5, Reg6, Reg7,Reg8)# Genero una Lista donde almaceno todas las regiones #--------------------------------------------------------------------------------------------------------------------- # BLOQUE 4: ANALISIS REGIONAL DE FRECUENCIAS Y ANALISIS DE TENDENCIA #--------------------------------------------------------------------------------------------------------------------- # PASO 8: DECLARACION DE VARIABLES PARA ALMACENAMIENTO DE RESULTADOS Regiones<-length(BaseRegiones) ResultadosSummaryStatistics<-array(0,dim=c(2400,7,Regiones))#2400=Cantidad de estaciones de la región que más estaciones tiene ResultadosSummaryStatisticsRegData<-array(0,dim=c(2400,7,Regiones))#2400=Cantidad de estaciones de la región que más estaciones tiene ResultadosRlmoments<-array(0,dim=c(5,Regiones))#5= L-momentos regionales ResultadosARFD<-array(0,dim=c(2400,Regiones))#2400=Cantidad estaciones máxima por región ResultadosARFH<-array(0,dim=c(3,Regiones))# 3= Indices de homogeneidad H1,H2,H3 ResultadosARFZ<-array(0,dim=c(5,Regiones))# 5= número de modelos de probabilidad a calcular su Bondad de Ajuste (glo, gev, gno, pe3, gpa) Resultadosrfitdist<-array(0,dim=c(1,Regiones))# 1=Un solo ajuste por región Resultadosrfitpara<-array(0,dim=c(5,Regiones))#5= cantidad de parametros de la Wakeby ResultadosRegionalQuantiles<-array(0,dim=c(19,Regiones))# 19=Máxima cantidad de cuantiles a calcular ResultadosRMAP<-array(0,dim=c(1,Regiones))# 1= Un solo valor de precipitacion media anual por region # PASO 9: ANALISIS REGIONAL DE FRECUENCIAS L-MOMENTOS PROPIAMENTE TAL, POR REGION for (z in 1:Regiones) { par(mfrow=c(1,2)) #Los comandos que siguen a continuacion permiten calcular todo el procedimiento paso a paso # Este primer comando calcula los l-momentos para las variables que van en las columnas designadas en Dataset[primera:última] # En caso de que se tengan más o menos, cambiar los valores. SummaryStatistics<-regsamlmu (BaseRegiones[[z]]) # El siguiente comando solamente convierte la salida anterior en un formato que puede ser leido por el comando siguiente. SummaryStatisticsRegData<-as.regdata(SummaryStatistics) # El siguiente comando grafica un Diagrama de L-momento ratios lmrd(SummaryStatisticsRegData) # Este comando calcula los L-momentos de la región constituida por las estaciones analizadas Rlmoments<-regavlmom(SummaryStatisticsRegData) # El siguiente comando agrega un punto rojo con los L-momentos regionales, al Diagrama de L-momentos-ratios lmrdpoints(Rlmoments, type="p", pch=22, col="red" ) # El siguiente comando calcula los estadísticos regionales, como el test de homogeneidad, las medidas de bondad de ajuste para varias distribuciones ARF<-regtst(SummaryStatisticsRegData, nsim=1000) # El siguiente comando ajusta el modelo de distribución seleccionado e identificado entre paréntesis. Usar el código del modelo apropiado. # Se declaran variables para almacenar los resultados de cada región, por separado # Se almacenan los resultados de Discordancia, Homegeneidad y Bondad de Ajuste a<-length(BaseRegiones[[z]]) #b<-nrow(BaseRegiones[[z]]) ResultadosSummaryStatistics[1:a,1:7,z]<-as.matrix(SummaryStatistics) #Se supone que aca almaceno todos los L-momentos ResultadosRlmoments[1:5,z]<-Rlmoments ResultadosARFD[1:a,z]<-ARF$D # Se almacenan las medidas de discordancia ResultadosARFH[1:3,z]<-ARF$H # Se almacenan las medidas de homogeneidad ResultadosARFZ[1:5,z]<-ARF$Z # Se almacenan las medidas de bondad de ajuste # PASO 10: SELECCION Y AJUSTE DEL MODELO DE DISTRIBUCION DE PROBABILIDAD rfit<-regfit(SummaryStatisticsRegData, "pe3") # Aca se selecciona Pearson III. Modificar según el modelo de mejor ajuste # El siguiente comando, calcula los cuantiles regionales para distintas probabilidad acumuladas RegionalQuantiles<-regquant(seq(0.05, 0.95, by=0.05), rfit) # Las siguientes tres lineas de comandos generan una gráfica de cuantiles rgc <- regqfunc(rfit)# Calcular Regional Growth Curve rgc(seq(0.05, 0.95, by=0.05)) curve(rgc, 0.01, 0.99, xlab="Non-exceedence Probability, F", ylab="Growth Curve") Resultadosrfitdist[z]<-rfit$dist # Se identifica la distribución utilizada Resultadosrfitpara[1:3,z]<-rfit$para # Se presentan los parametros de la distribución ajustada ResultadosRegionalQuantiles[1:19,z]<-RegionalQuantiles # Para cada región "z", almaceno sus resultados ResultadosRMAP[z]<-weighted.mean(SummaryStatisticsRegData[[3]],SummaryStatisticsRegData[[2]]) # Se calcula la precipitación media cada Región }# Cierra el ciclo for #cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #UNIR BASE DATOS REGIONAL TOTAL DE L-MOMENTOS CON BASE DATOS ESTACIONES. LUEGO EXPORTO Y PUEDO GENERAR LOS MAPAS DE L-MOMENTOS Total_Lmomentos<-as.data.frame(ResultadosSummaryStatistics[,,1]) colnames(Total_Lmomentos)[1]<-'id_estacion' colnames(Total_Lmomentos)[2]<-'n' colnames(Total_Lmomentos)[3]<-'Panalizada' colnames(Total_Lmomentos)[4]<-'L-CV' colnames(Total_Lmomentos)[5]<-'L-Skewness' colnames(Total_Lmomentos)[6]<-'L-Kurtosis' colnames(Total_Lmomentos)[7]<-'t5' BaseDatosTotal_Lmomentos<-merge(BaseConsolidada_sin_NA,Total_Lmomentos[,-2],by.x="id_estacion",by.y="id_estacion") # Acá uno las dos bases de datos a nivel de estaciones solamente. # ACA VOY A CALCULAR LA SIGNIFICANCIA DE LA PENDIENTE PARA EVALUAR SI HAY O NO TENDENCIA RegionTotalS<-sqldf("select id_estacion, anio, SumaLluviaAnual from BaseCompleta where LongRegistro>14") RegionTotalS_dat<-RegionTotalS[c("SumaLluviaAnual","anio")][,] RegionTotalS_fac<-factor(RegionTotalS["id_estacion"][,]) RegTotalS<-split(RegionTotalS_dat,RegionTotalS_fac)# estaciones<-length(RegTotalS) #.................................................................................... for (est in 1:estaciones){ Y<-RegTotalS[[est]]$SumaLluviaAnual/mean(RegTotalS[[est]]$SumaLluviaAnual) X<-RegTotalS[[est]]$anio-1900 fit = lm(Y ~ X) DWpvalue<-durbinWatsonTest(fit)$p lmp <- function (modelobject) { if (class(modelobject) != "lm") stop("Not an object of class 'lm' ") f <- summary(modelobject)$fstatistic p <- pf(f[1],f[2],f[3],lower.tail=F) s<-summary(modelobject)$coefficients[2] attributes(s) <- NULL attributes(p) <- NULL salida<-list(slope=s,Pvalue=p) return(salida) } RegTotalS[[est]]$Pendiente<-lmp(fit)$slope RegTotalS[[est]]$Pvalue<-lmp(fit)$Pvalue RegTotalS[[est]]$DW_Pvalue<-DWpvalue } RegTotalS2<-ldply(RegTotalS, data.frame) colnames(RegTotalS2)[1]<-'id_estacion' Pendiente_por_Estacion<-as.matrix(tapply(RegTotalS2[[4]],RegTotalS2[[1]],mean)) hist(Pendiente_por_Estacion) boxplot(Pendiente_por_Estacion) Pvalue_por_Estacion<-as.matrix(tapply(RegTotalS2[[5]],RegTotalS2[[1]],mean)) hist(Pvalue_por_Estacion) boxplot(Pvalue_por_Estacion) DW_por_Estacion<-as.matrix(tapply(RegTotalS2[[6]],RegTotalS2[[1]],mean)) hist(DW_por_Estacion) boxplot(DW_por_Estacion) id_estaciones<-levels(as.factor(RegTotalS2[,1])) BaseDatosTendenciaDW<-cbind(id_estaciones,Pendiente_por_Estacion,Pvalue_por_Estacion,DW_por_Estacion) colnames(BaseDatosTendenciaDW)[c(1,2,3,4)]<-c('id_estacion','Pendiente','Pvalue','D-W_Pvalue') BaseDatosTotal_Lmomentos<-merge(BaseDatosTotal_Lmomentos,BaseDatosTendenciaDW,by.x="id_estacion",by.y="id_estacion") # Acá uno las dos bases de datos a nivel de estaciones solamente. write.csv(BaseDatosTotal_Lmomentos, file = paste("BaseDatosTotal",ISO_PAIS,"_Lmomentos.csv",sep=''),row.names=FALSE) remove(BaseDatosTotal_Lmomentos) BaseDatosTotal_Lmomentos<-read.csv(paste("BaseDatosTotal",ISO_PAIS,"_Lmomentos.csv",sep='')) MAP<-BaseDatosTotal_Lmomentos[[9]] L_CV<-BaseDatosTotal_Lmomentos[[25]] L_Skewness<-BaseDatosTotal_Lmomentos[[26]] L_Kurtosis<-BaseDatosTotal_Lmomentos[[27]] plot(MAP,L_CV,main='L-CV versus MAP',xlab='Precipitación Media Anual [mm]',ylab='L-CV') plot(MAP,L_Skewness,main='L-Skewness versus MAP',xlab='Precipitación Media Anual [mm]',ylab='L-Skewness') plot(MAP,L_Kurtosis,main='L-Kurtosis versus MAP',xlab='Precipitación Media Anual [mm]',ylab='L-Kurtosis') #ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #------------------------------------------------------------------------------------------------------------------------ # BLOQUE 5: FUNCION DE AJUSTE DE CURVA L-MOMENTOS VS PRECIPITACION MEDIA ANUAL #------------------------------------------------------------------------------------------------------------------------ # PASO 11: DECLARACION DE VARIABLES PARA PROCESAR Y ALMACENAR COEFICIENTES DE AJUSTE DE CURVA RLCV <- ResultadosRlmoments[2,] RLSkewness<-ResultadosRlmoments[3,] RLKurtosis<-ResultadosRlmoments[4,] RMAP<-as.numeric(ResultadosRMAP) MAPvsLCV <- data.frame(RMAP,RLCV) MAPvsLSkewness<- data.frame(RMAP,RLSkewness) MAPvsLKurtosis<- data.frame(RMAP,RLKurtosis) # PASO 12: AJUSTE MEDIANTE ALTERNATIVAS DE OPTIMIZACION # ..................................................................................... # OPCION AJUSTE 1: Minimización mediante DEoptim. UTILIZAR ESTE METODO PMediaAnual<-RMAP LCVOBS<-RLCV LCVEST<-function(p) p[1]*exp(p[2]*PMediaAnual)+p[3] fun<-function(p) sum((LCVOBS-LCVEST(p))^2) ss <- DEoptim(fun, lower=c(0,-0.1,0), upper=c(0.3,0,0.2), control=list(trace=FALSE)) paLCV <- ss$optim$bestmem paLCV LSkOBS<-RLSkewness LSkEST<-function(p) p[1]*exp(p[2]*PMediaAnual)+p[3] fun<-function(p) sum((LSkOBS-LSkEST(p))^2) ss <- DEoptim(fun, lower=c(0,-0.1,0), upper=c(0.3,0,0.2), control=list(trace=FALSE)) paLSk <- ss$optim$bestmem paLSk LKurtOBS<-RLKurtosis LKurtEST<-function(p) p[1]*exp(p[2]*PMediaAnual)+p[3] fun<-function(p) sum((LKurtOBS-LKurtEST(p))^2) ss <- DEoptim(fun, lower=c(0,-0.1,0), upper=c(0.3,0,0.2), control=list(trace=FALSE)) paLKurt <- ss$optim$bestmem paLKurt #....................................................................................... # OPCION AJUSTE 2: Optimización mediante comando NLS (Non-linear Squares) #nlsfitLCV <- nls(RLCV~A*exp(B*RMAP)+C,data=MAPvsLCV, start=list(A=paLCV[1], B=paLCV[2], C=paLCV[3])) #nlsfitLSkewness <- nls(RLSkewness~A*exp(B*RMAP)+C,data=MAPvsLSkewness, start=list(A=paLSk[1], B=paLSk[2], C=paLSk[3])) #nlsfitLKurtosis <- nls(RLKurtosis~A*exp(B*RMAP)+C,data=MAPvsLKurtosis, start=list(A=paLKurt[1], B=paLKurt[2], C=paLKurt[3])) #pp<-seq(min(RMAP),max(RMAP),length=100) #plot(RMAP, RLCV, xlim=c(min(RMAP),max(RMAP)), ylim=c(min(RLCV),max(RLCV))) #lines(pp,predict(nlsfitLCV,list(RMAP=pp))) #plot(RMAP, RLSkewness, xlim=c(min(RMAP),max(RMAP)), ylim=c(min(RLSkewness),max(RLSkewness))) #lines(pp,predict(nlsfitLSkewness,list(RMAP=pp))) #plot(RMAP, RLKurtosis, xlim=c(min(RMAP),max(RMAP)), ylim=c(min(RLKurtosis),max(RLKurtosis))) #lines(pp,predict(nlsfitLKurtosis,list(RMAP=pp))) #summary(nlsfitLCV) #summary(nlsfitLSkewness) #summary(nlsfitLKurtosis) #.............................................................................................................. # OPCION AJUSTE 3: Minimización mediante comnado NLM (Non-Linear Minimization). Entrega los peores resultados. # Aca se presenta alternativa 2 para estimar mejor ajuste. #fnLCV <- function(p) sum((RLCV - p[1]*exp(p[2]*RMAP)+p[3])^2) #outLCV <- nlm(fnLCV, p = c(paLCV[1], paLCV[2], paLCV[3])) #outLCV$estimate #fnLSkewness <- function(p) sum((RLSkewness - p[1]*exp(p[2]*RMAP)+p[3])^2) #outLSkewness <- nlm(fnLSkewness, p = c(paLSk[1], paLSk[2],paLSk[3])) #outLSkewness$estimate #fnLKurtosis <- function(p) sum((RLKurtosis - p[1]*exp(p[2]*RMAP)+p[3])^2) #outLKurtosis <- nlm(fnLKurtosis, p = c(paLKurt[1], paLKurt[2], paLKurt[3])) #outLKurtosis$estimate #--------------------------------------------------------------------------------------------------------------------- # BLOQUE 6: GENERACION MAPA FRECUENCIAS Y PROBABILIDAD (PERIODO DE RETORNO) #--------------------------------------------------------------------------------------------------------------------- # PASO 13: IMPORTACION DE MAPA BASE TEMATICO DE LA VARIABILIDAD ESPACIAL DE LA VARIABLE PREDICTORA Mapa<-readGDAL("Mapa.tif") # Defino archivo de mapa base r<-raster(Mapa) # ESTE ES EL MAPA BASE, NO NECESARIAMENTE IGUAL AL MAPA DE VARIABLE FISICA ESPECIFICO projection(r) <- "+proj=latlong +ellps=WGS84" # Cambio proyección cartográfica # PASO 14: CALCULO DE LOS MAPAS DE L-MOMENTOS LCVmap<-paLCV[1]*exp(paLCV[2]*r)+paLCV[3] # Creo mapa de L-CV usando coeficientes de mejor ajuste LSmap<-paLSk[1]*exp(paLSk[2]*r)+paLSk[3] # Creo mapa de L-skewness usando coeficientes de mejor ajuste LKmap<-paLKurt[1]*exp(paLKurt[2]*r)+paLKurt[3] # Creo mapa de L-kurtosis usando coeficientes de mejor ajuste # PASO 15: CONVERSION DE FORMATO RASTER A MATRIZ PARA FACILIDAD DE CALCULO R<-as.matrix(r) # Este corresponde al MAPA DE VARIABLE FISICA ESPECIFICO. SOLO SE USARIA PARA CALCULO DE CUANTILES, NO DE PROBABILIDADES J<-as.matrix(LCVmap) K<-as.matrix(LSmap) L<-as.matrix(LKmap) # PASO 16: CALCULO DE LOS PARAMETROS DEL MODELO DE DISTRIBUCION. # Obtención de parámetros del(los) modelo(s) de probabilidad seleccionados(s). En este ejemplo, Pearson III. # Sólo usar una de las opciones y dejar el resto como comentario (#) Pearson3<-par.gamma((R/R),J,K) # Acá genero el mapa de parámetros para distribución Pearson según Viglione (alfa, beta,xi) (R/R es para crear un raster de 1s) #GenPar<-par.genpar((R/R),J,K) # Acá genero el mapa de parámetros para distribución Generalized Pareto según Viglione (alfa, beta,xi) (R/R es para crear un raster de 1s) #GEV<-par.GEV((R/R),J,K)# Acá genero el mapa de parámetros para distribución Generalized Extreme Value según Viglione (alfa, beta,xi) (R/R es para crear un raster de 1s) #LogNorm<-par.lognorm((R/R),J,K)# Acá genero el mapa de parámetros para distribución LogNormal según Viglione (alfa, beta,xi) (R/R es para crear un raster de 1s) #GenLogis<-par.genlogis((R/R),J,K)# Acá genero el mapa de parámetros para distribución Generalized Logistic según Viglione (alfa, beta,xi) (R/R es para crer un raster de 1s) #Kappa<-par.kappa((R/R),J,K,L)# Acá genero el mapa de parámetros para distribución Kappa según Viglione (alfa, beta,xi) (R/R es para crear un raster de 1s) # PASO 17: CALCULO DE LOS MAPAS DE FRECUENCIA/PERIODO DE RETORNO SEGUN CUANTIL #Acá, finalmente, obtengo el mapa de salida (probabilidad) o periodo de retorno (1/Dist) para un cuantil determinado (definido como un % de la VARIABLE FISICA ANALIZADA. SI ES LLUVIA EN JULIO, ES EL % DE ESA LLUVIA MEDIA DE JULIO) # Solo usar una de las opciones y dejar el resto como comentario (#) Cuantil<-0.4 # Como se esta analizando la precipitación anual, en este caso es 0.4 de la PMA. FreqMap<-F.gamma (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz #FreqMap<-F.genpar (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz #FreqMap<-F.GEV (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz #FreqMap<-F.lognorm (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz #FreqMap<-F.genlogis (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz #FreqMap<-F.kappa (Cuantil*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz ReturnMap<-(1/FreqMap)# Mapa de Periodo de Retorno formato Matriz MapaFrecuencia<-raster(FreqMap) MapaPeriodoRetorno<-raster(ReturnMap) # PASO 18: CALCULO DE LOS MAPAS DE FRECUENCIA SEGUN CANTIDAD DE LLUVIA #Acá, finalmente, obtengo el mapa de salida (probabilidad) o periodo de retorno (1/Dist) para un cuantil determinado (definido como un % de la VARIABLE FISICA ANALIZADA. SI ES LLUVIA EN JULIO, ES EL % DE ESA LLUVIA MEDIA DE JULIO) # Solo usar una de las opciones y dejar el resto como comentario (#) Lluvia<-500# Acá pregunto por la probabilidad de tener 200 mm FreqMap<-F.gamma ((Lluvia/R)*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz ReturnMap<-(1/FreqMap)# Mapa de Periodo de Retorno formato Matriz MapaFrecuencia<-raster(FreqMap) MapaPeriodoRetorno<-raster(ReturnMap) # PASO 19: CALCULO DE LOS MAPAS DE CANTIDAD DE LLUVIA SEGUN PROBABILIDAD #Acá, finalmente, obtengo el mapa de salida (probabilidad) o periodo de retorno (1/Dist) para un cuantil determinado (definido como un % de la VARIABLE FISICA ANALIZADA. SI ES LLUVIA EN JULIO, ES EL % DE ESA LLUVIA MEDIA DE JULIO) # Solo usar una de las opciones y dejar el resto como comentario (#) Fx<-0.2# Acá pregunto por la probabilidad de tener 200 mm CantidadLLuvia<-invF.gamma (Fx*(R/R), Pearson3$xi, Pearson3$beta, Pearson3$alfa)# Mapa de Probabilidad formato Matriz MapaLluvia<-raster(CantidadLluvia) #PASO 18: APLICACION DE FORMATO RASTER A LOS MAPAS DE FRECUENCIA/PERIODO DE RETORNO projection(MapaFrecuencia) <- "+proj=latlong +ellps=WGS84" # Cambio proyección cartográfica #dim(MapaFrecuencia)<-(973,456) # Asigno al mapa de salida las mismas dimensiones que el mapa base xmin(MapaFrecuencia)<-xmin(r) # Asigno al mapa de salida los mismos límites que el mapa base xmax(MapaFrecuencia)<-xmax(r) # Asigno al mapa de salida los mismos límites que el mapa base ymin(MapaFrecuencia)<-ymin(r) # Asigno al mapa de salida los mismos límites que el mapa base ymax(MapaFrecuencia)<-ymax(r) # Asigno al mapa de salida los mismos límites que el mapa base projection(MapaPeriodoRetorno) <- "+proj=latlong +ellps=WGS84" # Cambio proyección cartográfica #dim(MapaPeriodoRetorno)<-(973,456) # Asigno al mapa de salida las mismas dimensiones que el mapa base xmin(MapaPeriodoRetorno)<-xmin(r) # Asigno al mapa de salida los mismos límites que el mapa base xmax(MapaPeriodoRetorno)<-xmax(r) # Asigno al mapa de salida los mismos límites que el mapa base ymin(MapaPeriodoRetorno)<-ymin(r) # Asigno al mapa de salida los mismos límites que el mapa base ymax(MapaPeriodoRetorno)<-ymax(r) # Asigno al mapa de salida los mismos límites que el mapa base projection(MapaLluvia) <- "+proj=latlong +ellps=WGS84" # Cambio proyección cartográfica #dim(MapaPeriodoRetorno)<-(973,456) # Asigno al mapa de salida las mismas dimensiones que el mapa base xmin(MapaLluvia)<-xmin(r) # Asigno al mapa de salida los mismos límites que el mapa base xmax(MapaLluvia)<-xmax(r) # Asigno al mapa de salida los mismos límites que el mapa base ymin(MapaLluvia)<-ymin(r) # Asigno al mapa de salida los mismos límites que el mapa base ymax(MapaLluvia)<-ymax(r) # Asigno al mapa de salida los mismos límites que el mapa base # PASO 19: PLOTEO DE LOS MAPAS DE FRECUENCIA/PERIODO DE RETORNO # Ploteo plot(r) # Aca ploteo mapa base plot(LCVmap) # Ploteo mapa de L-Cv plot(LSmap) # Plotep mapa de L-Skewness plot(LKmap) # Ploteo mapa de L-Kurtosis plot(MapaFrecuencia) plot(MapaPeriodoRetorno) # PASO 20: EXPORTACION DE LOS MAPAS A FORMATO GEOTIFF # Exportación de mapas writeRaster(MapaFrecuencia,filename="MapaFrecuencia.tif",format="GTiff",overwrite=TRUE,NAflag=-999) writeRaster(MapaPeriodoRetorno,filename="MapaPeriodoRetorno.tif",format="GTiff", overwrite=TRUE,NAflag=-999) writeRaster(MapaLluvia,filename="MapaLluvia.tif",format="GTiff", overwrite=TRUE,NAflag=-999)