#ggplot2 #R Comparación de modelos de pronósticos de CO y NOx en Buenos Aires, Bogotá y Santiago, con regresiones y redes neuronales en R, ggplot2 y openair

-
Para el congreso Argentina Ambiental 2012 he estado terminando un resumen extendido, y me tiene muy contento el hecho de aplicar estadística dentro del ámbito de los sistemas ambientales.

Para este congreso realicé un pronóstico de CO y  NOx para las ciudades de Buenos Aires, Bogotá y Santiago, comparando los resultados de una regresión y redes neuronales. A continuación presentaré dos gráficos que comparan los resultados de forma muy simple y didáctica.

Entonces, con las técnicas que no voy a presentar acá genere las regresiones y redes neuronales, y obtuve un CO y NOx pronosticado con cada técnica para cada ciudad. Pensé en graficar con ggplot2 aplicando un factor específico, y se me ocurrió en excel clasificar las  el valor absoluto de la diferencia de los pronósticado con lo original, según percentiles del original, así:

abs(diferencia)<P(20)=1; es decir, el valor absoluto de la diferencia debe cae dentro del percentil 20 del valor original para que de el valor de 1, y entonces generé la escala

abs(diferencia)

  • <P(20)=1
  • P(20)-P(40)=0.8
  • P(40)-P(60)=0.6
  • P(60)-P(80)=0.4
  • >P(80)=0.2

Y la fórmula en Excel 2010 es la siguiente:

SI(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.2);1;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.2);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.4)));0.8;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.4);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.6)));0.6;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.6);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.8)));0.4;SI(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.8);0.2;0)))))

y ahí tenemos nuestro factor


Luego:


setwd("C:/r")
library(openair)
#Uso la librería openair ya que es ideal para trabajar con datos en formato típico
#attach es para manejar directamente las variables, pero acá no es necesario
aire1<-import("rpud-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire1)
names(aire1)
aire2<-import("rpud-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire2)
names(aire2)
aire3<-import("rbai-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire3)
names(aire3)
aire4<-import("rbai-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire4)
names(aire4)
aire5<-import("rbog-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire5)
names(aire5)
aire6<-import("rbog-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire6)
names(aire6)
library(ggplot2)
Y ahora podemos graficar






> p7 <-ggplot(aire6)+aes(x=bog_nox_regpro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_reg),colour = factor(bog_dif_co_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p10 <- ggplot(aire6)+aes(x=bog_nox_neupro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_neu),colour = factor(bog_dif_co_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p8 <-ggplot(aire4)+aes(x=bai_nox_reg_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_reg),colour = factor(bai_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 140,y = 300, label = "R2 = 0.40", colour = "black", size=8)+ opts(title="NOx Buenos Aires Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p11 <-ggplot(aire4)+aes(x=bai_nox_neu_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_neu),colour = factor(bai_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 100,y = 300, label = "R2 = 0.33", colour = "black", size=8)+ opts(title="NOx Buenos Aires Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p9 <-ggplot(aire2)+aes(x=p_nox_reg_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_reg),colour = factor(p_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.73", colour = "black", size=8)+ opts(title="NOx Santiago Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p12 <- ggplot(aire2)+aes(x=p_nox_neu_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_neu),colour = factor(p_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.80", colour = "black", size=8)+ opts(title="NOx Santiago Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)"
> p7 <-ggplot(aire6)+aes(x=bog_nox_regpro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_reg),colour = factor(bog_dif_co_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.46", colour = "black", size=8)+ opts(title="NOx Bogotá Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p10 <- ggplot(aire6)+aes(x=bog_nox_neupro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_neu),colour = factor(bog_dif_co_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p8 <-ggplot(aire4)+aes(x=bai_nox_reg_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_reg),colour = factor(bai_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 140,y = 300, label = "R2 = 0.40", colour = "black", size=8)+ opts(title="NOx Buenos Aires Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p11 <-ggplot(aire4)+aes(x=bai_nox_neu_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_neu),colour = factor(bai_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 100,y = 300, label = "R2 = 0.33", colour = "black", size=8)+ opts(title="NOx Buenos Aires Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p9 <-ggplot(aire2)+aes(x=p_nox_reg_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_reg),colour = factor(p_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.73", colour = "black", size=8)+ opts(title="NOx Santiago Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p12 <- ggplot(aire2)+aes(x=p_nox_neu_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_neu),colour = factor(p_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.80", colour = "black", size=8)+ opts(title="NOx Santiago Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")


Y con esta fórmula sensacional de la receta de cocina para R, el Cookbook for R


# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

obtenemos

multiplot(p7, p8, p9, p10, p11, p12,cols=3)




multiplot(p7, p8, p9, p10, p11, p12,cols=3)




Referencias:
R Development Core Team (2011). R: A language and environment for statistical computing. R
  Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL
  http://www.R-project.org/.
H. Wickham. ggplot2: elegant graphics for data analysis. Springer New York, 2009.


Teetor Paul. Cookbook for R. Pub. O'Really Media. 2011. 978-0-596-80915-7

Comentarios

Publicar un comentario

Entradas populares de este blog

ERROR: Grib2 file or date problem, stopping in edition_num. Can't ungrib.exe SOLVED

-
Imagen
Well, it seems that this problem can be very common. Usually, you have a directory for downloading the meteorology from NCEP/UCAR. Probably, you will have the script for downloading the data in the same directory. If you link the grib files based on the name there won't be any problems ./link_grib.csh /path/GFS* However, it could be possible to link the whole directory. You might do this when you have a whole year of meteorological data. ./link_grib.csh /path/ If you have scripts for downloading the met data in the same directory, THEN you will have the message: ERROR: Grib2 file or date problem, stopping in edition_num To avoid this error, just move the script for downloading the met data elsewhere and link the directory again. Finally, do ungrib.exe UPDATE 2018-12/15 (I mean 15 of December of 2018) I had the same problem again, and with the same data! and this solution did not work till i did this: ./link_grib.csh /path/GFS Yes only the first letters, in
Leer más

ARCGIS: Cálculo de pendiente entre puntos en el espacio xyz

-
Imagen
Se tienen puntos en el espacio con coordenadas xyz. Estas corresponden a coordenadas UTM huso 19, en metros. Pendiente = Distancia Vertical / Distancia Horizontal  Distancia Vertical = z2 - z1 Distancia Horizontal = d = [(y2-y1)^2 + (x2-x1)^2]^0.5 Pendiente = (z2 - z1) / [(y2-y1)^2 + (x2-x1)^2]^0.5 Pendiente porcentual = Pendiente *100 Entonces; x1 y1 z1 x2 y2 z2 pendiente% 344150.34 6302722 520.897 344911.2 6302394 530.5 1.159 344150.34 6302722 520.897 343774.7 6302378 521.858 0.188 345811.8 6302371 536 345893.9 6302547 533 -1.544 Sirve para calcular pendiente de caminos, ríos, etc :) En este caso lo hice para una modelación de flujo vehicular de Santiago Estraus , por lo que la diferencias de las alturas está según el sentido del flujo, desde la altura del nodo inicial al nodo final, por eso hacia el sector oriente se aprecia mayor pendiente. Además, se aprecian
Leer más

WRF open_aux_u : error opening auxinput5_d02_2013-12-31_00:00:00 for reading. 100

-
Imagen
So, i got this error: open_aux_u Timing for Writing /scr4/sibarra/masp/health/wrfo for domain        2:    4.14512 elapsed seconds open_aux_u : error opening auxinput5_d02_2013-12-31_00:00:00 for reading. 100 pretty weird because it a met run. According to http://forum.wrfforum.com/viewtopic.php?f=47&t=10003, they recommend diminishing the number of processors. I ran with 32 and got stuck I ran with 16 and appear this error open_aux_u : error opening auxinput5_d02_2013-12-31_00:00:00 for reading. 100 I'm running with  8 and still running. 
Leer más