04-capitulo3.Rmd

# Asimilacion de radianzas del satelite geoestacionario GOES-16 en cielos despejados {#cap-5-goes}

En este capítulo se evalúa el impacto que genera la asimilación de observaciones en cielos despejados del sensor ABI a bordo del satélite GOES-16 tanto en el análisis como en pronósticos determinísticos de alta resolución para el caso de estudio del 22 de noviembre de 2018. Estas observaciones son de particular interés en la región de Sudamérica por su alta resolución espacial y frecuencia temporal. Usaremos las observaciones de los canales sensibles al vapor de agua que brindan información en 3 niveles de la atmósfera, algo que no era posible con los sensores previos a bordo de los satélites GOES. Esperamos que la asimilación de estas observaciones contribuya a mejorar la representación de los procesos convectivos y de mesoescala en pronósticos a corto plazo. En particular, los pronósticos en alta resolución permitiran evaluar el desarrollo de la convección que produce el modelo de manera explícita.

Por ser un satélite y sensor nuevos, la última versión de GSI (V3.7) que fue publicada en julio de 2018 no incluye la posibilidad de asimilar estas observaciones. Fue necesario, entonces, modificar el código fuente del sistema de asimilación para incluir las rutinas necesarias utilizando como base los trabajos previos de @lee2019 y @liu2019. También por esta razón, el capítulo incluirá el análisis de experimentos para definir la mejor configuración posible y evaluar técnicamente la asimilación de observaciones de ABI.

Trabajos previos evaluaron el impacto de asimilar observaciones de los canales asociados al vapor de agua de ABI en escala regional. En líneas general estos trabajos (@lee2019, @jones2020, entre otros) utilizan las observaciones de los 3 canales en conjunto. Sin embargo @honda2018, utilizaron observaciones del sensor Advanced Himawari Imager (AHI), de iguales características que ABI, y encontraron mejores resultados al asimilar los canales 8 (6.2 $\mu m$) y 10 (7.3 $\mu m$). También observaron que el canal 9 (6.9 $\mu m$) está muy correlacionado con los dos primeros y por lo tanto no aporta información independiente en la asimilación. Por otro lado @lee2019 observaron que el impacto era mayor al asimilar observaciones de los 3 canales sobre el continente y el mar al mismo tiempo, incorporando así más observaciones en el análisis. En el mismo sentido, observaron que el incremento del análisis disminuía considerablemente al remover el canal 10 y en algunos casos también el canal 9, lo que indica que estos canales producen un impacto importante.

En la asimilación de datos y en particular de radianzas de satélites, muchas observaciones con alta resolución espacial deben ser descartadas porque podrían tener un impacto negativo en el análisis al no ser independientes entre si [@lazarus2010; @janjic2018]. En los esquemas de asimilación actuales se asume que la matriz de covarianza de los errores de las observaciones $R$ es diagonal, es decir, no hay correlación entre los errores de las distintas observaciones y cada observación es independiente. Sin embargo, las observaciones de satélites pueden estar muy correlacionadas debido a su resolución espectral y espacial. En este contexto se utilizan estrategias como el thinning para disminuir la resolución de las observaciones asumiendo que la correlación decae con la distancia y conservando la información que aportan. Por ejemplo @jones2020 aplican la técnica de thinning para llevar las observaciones con resolución de 2 km a una resolución de 15 km y así evitar el impacto de la correlación espacial entre las observaciones.  

## Metodología

### Asimilación de observaciones de GOES-16 {#asim-abi}

La asimilación directa de radianzas de ABI requiere las mismas consideraciones descriptas previamente en la Sección \@ref(asim-rad) en cuanto a control de calidad, corrección de bias, thinning e identificación de pixeles nubosos. En esta Sección se describen las decisiones metodológicas y técnicas que se tomaron a la hora de asimilar las observaciones, en base a las pruebas preliminares realizadas. 

La asimilación de radianzas en *all sky*, es decir tanto para cielos claros como nubosos es todavía un importante desafío ya que estos dos tipos de observaciones requieren un control de calidad y estimaciones de errores específicos. En esta primera aproximación a la asimilación de observaciones de ABI solo se utilizaron observaciones en cielos despejados. Para detectar y filtrar los pixeles nubosos utilizamos la mascara de nubes ACM [por *ABI Cloud Mask,* @heidinger2013] disponible con la misma frecuencia que las observaciones. Esta mascara de nubes da información binaria (cielo despejado o cielo nuboso) y tienen una exactitud del 87% [@heidinger2013]. La clasificación es realizada utilizando información de 10 canales de ABI e incluye criterios espectrales, espaciales y temporales. En las Figuras \@ref(fig:cld-msk) se presentan 2 ejemplos de máscara de nubes a las 00 UTC del 20 de noviembre y a las 18 UTC del 22 de noviembre. 

(ref:cld-msk) Ubicación de píxeles nubosos (gris) según la máscara de nubes generada para GOES-16, durante el a) 06 UTC del 20 de noviembre y b) 18 UTC del 22 de noviembre. 

```{r cld-msk, fig.cap="(ref:cld-msk)", fig.width=6, fig.height=3}
#dominio de interés
xlim <- c(2600, 3950)
ylim <- c(3700, 4750)

files_acm <- Sys.glob(here("data/derived_data/observations/OR_ABI-L2-ACMF-M3_G16_*"))


cld_msk <- map(files_acm, function(f){
  
  nc <- ncdf4::nc_open(f)
  
  ReadNetCDF(f, vars = c('BCM', 't'), 
             subset = list(x = xlim,  y = ylim)) %>% 
    .[, t := as_datetime(t, origin = "2000-01-01 12:00:00", tz = "UTC")] %>% 
    .[, c("lon", "lat") := goes_projection(x, y, nc)]
  
}) %>% 
  rbindlist()


cld_msk %>% 
  .[, t := factor(t)] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = factor(BCM))) +
  scale_color_manual(values = c("white", "darkgrey"), guide = NULL) +
  geom_mapa() +
  facet_wrap(~t, labeller = labeller(t = c("2018-11-20 00:05:52" = "6 UTC, 20 de Nov",
                                           "2018-11-22 18:05:54" = "18 UTC, 22 de Nov"))) +
  tagger::tag_facets() +
  labs(x = NULL, y = NULL) +
  theme_minimal(base_size = 9)

```

Las radianzas de satélite pueden tener errores sistemáticos asociados a las condiciones de la atmósfera, a la geometría de la observación y condiciones específicas de los sensores que deben ser corregidos previo a la asimilación. Por este motivo se generó un primer experimento de asimilación para evaluar la magnitud de estos errores. Esta primera prueba utiliza la misma configuración de modelo y ensamble que los experimentos discutidos previamente e incluye la asimilación de los 3 canales de ABI sensibles al vapor de agua en conjunto.

La Figura \@ref(fig:oma-dist) muestra la distribución de la diferencia entre las observaciones y el campo preliminar (OmF) y entre las observaciones y el análisis (OmA) para todo el periodo. Si las observaciones no tienen bias la distribución OmF estará centrada en cero. En este caso solo el canal 9 presenta un bias positivo pequeño del orden de 0.3 K (Tabla \@ref(tab:oma-tabla)). Por esta razón y siguiendo trabajos previos [por ejemplo @honda2018] en los experimentos discutidos en este capítulo no se aplica una corrección del bias a estas observaciones. 

La Figura \@ref(fig:oma-dist) y la Tabla \@ref(tab:oma-tabla) muestran además otro resultado importante. Si la asimilación de las observaciones se realizó correctamente es esperable que la diferencia OmA esté más cerca de cero que OmF y que su distribución tenga menos varianza. Esto se observa para los 3 canales asimilados y confirma que la implementación técnica dentro del sistema funciona correctamente. Notar que los cambios abruptos y simétricos en las curvas de distribución se deben al chequeo preliminar que se realiza como parte de los controles de calidad y que elimina aquellas observaciones que tienen un diferencia mayor a 2.2 K con respecto del campo preliminar. 

(ref:oma-dist) Distribución de la diferencia entre las observaciones y el campo preliminar (OmF) y la observación y el análisis (OmA) para los 3 canales de ABI y a lo largo de todo el periodo de prueba. 

```{r oma-dist, fig.cap="(ref:oma-dist)", fig.align='center'}
# path_data <- "/home/paola.corrales/datosmunin3/EXP/"
# 
# files_gue <- Sys.glob(paste0(path_data, "E10/ANA/*/diagfiles/asim_abi_g16_*.ensmean"))
# files_ana <- Sys.glob(paste0(path_data, "E10/ANA/*/diagfiles_ana/asim_abi_g16_*.ensmean"))
# 
# d <- map2(files_gue, files_ana, function(fg, fa) {
#   
#   meta <- unglue::unglue(fg, "/home/paola.corrales/datosmunin3/EXP/E10/ANA/{date}/diagfiles/asim_abi_g16_{date2}.ensmean")
#   gue <- read_diag_rad(fg, "E10") %>% 
#     .[, channel := channel + 6] %>% 
#     .[channel %in% c(8:10) & qc == 0, .(channel, lat, lon, tb_obs, tbc)] %>% 
#     .[, run := "guess"]
#   
#   ana <- read_diag_rad(fa, "E10") %>% 
#     .[, channel := channel + 6] %>% 
#     .[channel %in% c(8:10) & qc == 0, .(channel, lat, lon, tb_obs, tbc)] %>% 
#     .[, run := "ana"]
#   
#   rbind(ana, gue) %>% 
#     .[, fecha := meta[[1]][["date"]]] %>% 
#     # .[, ":="(lon.box = cut_round(lon, breaks = seq(284, 309, 2.5)),
#              # lat.box = cut_round(lat, breaks = seq(-42, -17, 2.5)))] %>% 
#     .[]
#   
# }) %>% 
#   rbindlist()
# write_csv(d, here("data/derived_data/sample_obs/E10_omb_oma.csv"))

diag_E10 <- fread(here("data/derived_data/sample_obs/E10_omb_oma.csv"))

diag_E10 %>% 
  ggplot(aes(tbc)) +
  geom_density(aes(color = run), key_glyph = draw_key_path) +
  scale_x_continuous(limits = c(-5, 5)) +
  scale_color_manual(values = c("#1f78b4", "#33a02c"), labels = c("OmA", "OmF")) +
  labs(x = "OmF / OmA (K)",
       y = "Densidad de probabilidad", 
       color = NULL) +
  facet_wrap(~channel, labeller = labeller(channel = c("8" = "Canal 8",
                                                       "9" = "Canal 9",
                                                       "10" = "Canal 10"))) +
  tagger::tag_facets() +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

```{r oma-tabla}
diag_E10[, .(rms = mean(tbc)), by = .(channel, run)] %>% 
  dcast(run ~ channel) %>% 
  .[, run := c("OmA", "OmF")] %>% 
  arrange(desc(run)) %>% 
  kbl(booktabs = TRUE, col.names = c("", "Canal 8", "Canal 9", "Canal 10"),
      digits = 4, 
      align = "lrrr",
      caption = "Observación menos análisis o campo preliminar medio calculado sobre todo el periodo en K.") %>% 
  kable_classic_2(full_width = TRUE) 
```


Debido al desarrollo de actividad convectiva en el dominio durante el periodo simulado, la cantidad de observaciones en cielos despejados varía considerablemente (Figuras \@ref(fig:cld-msk)a y b). La Figura \@ref(fig:obs-rad)a muestra el número de radianzas asimiladas en cada ciclo de asimilación donde se aprecia una disminución del número de observaciones en las últimas 24 hs del periodo debido a la presencia de nubes. Sin embargo la cantidad de observaciones que aporta ABI es casi siempre entre 1 y 2 órdenes de magnitud mayor que el número de observaciones de satélites polares. En este caso ABI utiliza un thinning de 10 km (ver Sección \@ref(thinning)) pero aún cuando se utiliza un thinning de 30 km, el número de observaciones de ABI asimiladas es siempre mayor a lo que pueden aportar los satélites polares. Las observaciones de los canales 8, 9 y 10 de ABI, además, se distribuyen de manera homogénea en el dominio y cubren gran parte de los niveles verticales de la atmósfera. La Figura \@ref(fig:obs-rad)c muestra la distribución promedio de observaciones de cada canal con la altura. Si bien el nivel donde se observa el máximo de observaciones para cada canal es similar a lo esperado según sus funciones de peso (Figura \@ref(fig:wf-goes)), las observaciones se encuentran distribuidas en un rango amplio de alturas, ya que la posición del máximo de la función de peso en cada caso depende de las condiciones particulares de la atmósfera en cada punto que difiere de la atmósfera estándar utilizada para el cálculo de las funciones de peso teóricas.

(ref:obs-rad) a) Número de radianzas asimiladas en cada ciclo, b) número medio de radianzas asimiladas por ciclo en capas verticales de 50 hPa de espesor para los experimentos RAD (radianzas de satélites polares, círculos rojos) y ABI (radianzas de del sensor ABI, cuadrados turquesas) y c) número medio de radianzas asimiladas por ciclo en capaz verticales de 50 hPa de espesor para los 3 canales de ABI.

```{r obs-rad, fig.cap="(ref:obs-rad)", fig.height=4, fig.width=6}

files <- Sys.glob(here("data/derived_data/omb_diagfiles/E10/asim*ensmean"))

files <- files[!str_detect(files, "conv")]

rad <- map(files, function(f){ 
  
  rad <- read_diag_rad(f, "E10") %>% 
    .[, channel := fifelse(sensor == "abi_g16", channel + 6, channel)] %>% 
    satinfo[., on = c("sensor", "channel")] 
  
  rad <- rad[, .(sensor, peakwt, iuse, error, exp, date, qc, channel)] %>% 
    setnames(c("peakwt", "iuse", "error"), c("pressure", "usage.flag", "rerr"))
  
  rad_date <- rad %>% 
    # rbind(., rad) %>% 
    .[usage.flag >= 1 & rerr != 1.0e+10 & qc == 0] %>% 
    .[, .(count_obs = .N), by = .(sensor, exp, date)] 
  
  rad_lev <- rad %>% 
    .[usage.flag >= 1 & rerr != 1.0e+10 & qc == 0] %>% 
    .[, ":="(lev.box = cut_round(pressure, breaks = seq(0, 1050, 50)))] %>% 
    .[, .(count_obs = .N), by = .(sensor, channel, exp, lev.box, date)]
  list(rad_date = rad_date, rad_lev = rad_lev)
  
}) %>% 
  reduce(function(x, y) list(rad_date = rbind(x$rad_date, y$rad_date),
                             rad_lev = rbind(x$rad_lev, y$rad_lev)))


rad$rad_date %>% 
  .[, exp := fifelse(str_detect(sensor, "abi"), "ABI", "RAD")] %>% 
  .[, .(count_obs = sum(count_obs)), by = .(exp, date)] %>%
  ggplot(aes(date, count_obs)) +
  geom_line(aes(color = exp), size = 0.2) +
  geom_point(aes(color = exp, shape = exp), size = 1) +
  scale_color_manual(values = c("#047994", "#CC3311")) +
  scale_shape_manual(values = c(15, 16)) +
  scale_y_log10() +
  scale_date() +
  labs(y = "Número de obs por ciclo", x = "Hora (UTC)",
       color = NULL, shape = NULL) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  
  rad$rad_lev %>% 
  .[, exp := fifelse(str_detect(sensor, "abi_g16"), "ABI", "RAD")] %>% 
  .[, .(obs_cycle = mean(count_obs, na.rm = TRUE)), by = .(exp, lev.box)] %>% 
  ggplot(aes(lev.box, obs_cycle)) +
  geom_line(aes(color = exp), size = 0.2) +
  geom_point(aes(color = exp, shape = exp), size = 1) +
  scale_color_manual(values = c("#047994", "#CC3311")) +
  scale_shape_manual(values = c(15, 16)) +
  scale_y_log10() +
  scale_x_reverse(limits = c(1050, 50)) +
  coord_flip() +
  labs(color = NULL, shape = NULL,
       x = "Nivel de presión (hPa)", y = "Promedio de \nobs por ciclo") +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  
   rad$rad_lev %>% 
  .[, exp := fifelse(str_detect(sensor, "abi_g16"), "ABI", "RAD")] %>% 
  .[ exp == "ABI"] %>% 
  .[, .(obs_cycle = mean(count_obs, na.rm = TRUE)), by = .(exp, channel, lev.box)] %>% 
  ggplot(aes(lev.box, obs_cycle)) +
  ggnewscale::new_scale_color() +
  geom_line(aes(color = factor(channel)), size = 0.2) +
  geom_point(aes(color = factor(channel), shape = factor(channel)), size = 1) +
  # scale_color_manual(values = c("#047994", "#CC3311")) +
  scale_color_viridis_d(labels = c("Canal 8", "Canal 9", "Canal 10")) +
  scale_shape_manual(values = c(15, 16, 17), labels = c("Canal 8", "Canal 9", "Canal 10")) +
  scale_y_log10() +
  scale_x_reverse(limits = c(1050, 50)) +
  coord_flip() +
  labs(color = NULL, shape = NULL,
       x = NULL, y = "Promedio de \nobs por ciclo") +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  
  plot_layout(ncol = 3, widths = c(4, 1.25, 1.2), guides = "collect") & 
  plot_annotation(tag_levels = 'a', tag_suffix = ")") &
  theme(plot.tag = element_text(size = 8), 
        legend.position = "bottom",
        legend.spacing = unit(c(0, 0, 0, 0), "cm"),
        panel.background = element_rect(fill = "#fbfbfb", color = NA),
        plot.margin = unit(c(0, 0, 0.5, 0), "cm"))

```

### Configuracion de los experimentos

En este capítulo se analizan una serie de experimentos, primero para evaluar la sensibilidad a distintas combinaciones de canales y a la resolución del thinning y finalmente para comparar el impacto de la asimilación de observaciones de ABI en conjunto o no con observaciones de satélites polares. Los experimentos SATWND y RAD fueron analizados previamente en el Capítulo \@ref(cap-3-analisis). SATWND asimila observaciones convencionales, EMA y vientos derivados de satélite y RAD suma a lo anterior radianzas en aire claro de sensores a bordo de satélites polares. RAD+ABI incluye además observaciones de los 3 canales de vapor de agua de ABI con un thinning de 10 km. Finalmente el conjunto de experimentos ABI, descriptos en la Tabla \@ref(tab:exp-rad), incluyen radianzas de ABI pero no de satélites polares, y usan distintas combinaciones de canales y resolución de thinning para evaluar la sensibilidad de la asimilación de distintos canales y la resolución del thinning. 


```{r exp-rad}
tribble(
  ~exp, ~obs, ~canales, ~thinning,
  "SATWND", "--", "--", "--", 
  "RAD", "Radianzas de satélites polares", "--", "--",
  "RAD+ABI", "Radianzas de satélites polares y ABI", "8, 9, 10 ", "10 km ", 
  "ABI_ch890_th10 (ABI)", "Radianzas de ABI", "8, 9, 10 ", "10 km ", 
  "ABI_ch890_th30", "Radianzas de ABI", "8, 9, 10 ", "30 km ",
  "ABI_ch90_th30", "Radianzas de ABI", "9, 10 ", "30 km ",
  "ABI_ch0_th30", "Radianzas de ABI", "10 ", "30 km ",
) %>% 
  kbl(booktabs = TRUE, col.names = c("Experimento", "Radianzas asimiladas", "Canales ABI", "Thining ABI"),
      align = "llcc",
      caption = "Experimentos realizados. Todos los experimentos incluyen observaciones convencionales, EMA y vientos derivados de satélite. ") %>% 
  kable_classic_2(full_width = TRUE) %>% 
  kable_styling(font_size = 9) %>% 
  column_spec(1, width = "11em") %>%
  column_spec(2, width = "16em") %>%
  #   column_spec(2:3, width = "2.5em", latex_valign = "m") %>% 
  column_spec(3:4, width = "4em", latex_valign = "m")

```


Todos los experimentos previamente descriptos usan la configuración del ensamble multifísica y del sistema de asimilación explicada en la Sección \@ref(configmodelo). De la misma manera los ciclos de análisis horarios se realizan desde las 18 UTC del 20 de noviembre hasta las 12 UTC del 23 de noviembre para cubrir el desarrollo y maduración del SCM. 

Para estudiar el impacto de la asimilación de radianzas también se generaron pronósticos determinísticos en alta resolución inicializados a las 00 del 22 de noviembre (Figura \@ref(fig:cycle-fcst)). Debido al altísimo costo computación que requiere una simulación en alta resolución, no fue posible generar pronósticos por ensambles como en el Capítulo \@ref(cap-4-pronosticos) y entre otras cosas analizar los pronósticos de precipitación desde el punto de vista probabilístico. Sin embargo, estos pronósticos permiten analizar el desarrollo convectivo que el modelo resuelve explicitamente en mayor detalle, característica que puede mejorar con la asimilación de observaciones de ABI. En la Figura \@ref(fig:dominio-det) se muestra el dominio *d01* con 10 km de resolución (línea negra), idéntico al dominio utilizado para el resto de las simulaciones numéricas y el dominio anidado *d02* con 2 km de resolución (línea turquesa). Las condiciones iniciales están dadas por los análisis de los distintos experimentos y las condiciones de borde son generadas a partir del GFSF determinístico. Además las variables atmosféricas para el *d02* fueron inicializadas haciendo un escalado de la información del análisis. Todos los pronósticos utilizan las mismas parametrizaciones: del modelo de superficie terrestre [Noah-MP, @chen2001], de microfísica [esquema de un solo momento de 6 clases del WRF, @hong2006a] de procesos radiativos [esquema de onda corta y onda larga del RRTMG, @iacono2008], YSU [@hong2006] para capa límite y KF [@kain2004] para covección (solo aplicado al *d01*).

Además de los pronósticos inicializados a partir de los análisis, se generó un pronóstico inicializado a partir del análisis global de GFS, es decir, sin asimilación de datos regional. Este pronóstico tiene la misma configuración y usa el mismo conjunto de parametrizaciones que los pronósticos descriptos previamente. 

(ref:dominio-det) Dominio de baja resolución (d01, 10 km) marcado con un recuadro negro y dominio anidado de alta resolución (d02, 2 km) marcado en color turquesa. 

```{r dominio-det, fig.cap="(ref:dominio-det)", fig.width=3, fig.height=3, fig.align="center"}
# ReadNetCDF("/home/paola.corrales/datosmunin3/EXP/E10/FCST/2018112200_det/wrfout_d02_2018-11-22_00:00:00", vars = c("XLAT", "XLONG")) %>% 
#   setnames(c("XLAT", "XLONG"), c("lat", "lon")) %>% 
#   write_csv("/home/paola.corrales/datosmunin3/EXP/paper_ana_20181122-derived_data/sample_obs/dominio_d02.csv")

square_d02 <- fread(here("data/derived_data/sample_obs/dominio_d02.csv")) %>% 
  .[, square1 := west_east %in% range(west_east) , by = .(south_north)] %>% 
  .[, square2 := south_north %in% range(south_north) , by = .(west_east)] %>% 
  .[square1 | square2] 


dominio <- fread(here("data/derived_data/sample_obs/dominio_hgt.csv")) %>% 
  .[, c("x", "y") := wrf_project(lon, lat)]

dominio %>% 
  ggplot(aes(x, y)) +
  geom_contour_fill(aes(z = hgt), proj = norargentina_lambert,
                    breaks = seq(0, 6000, 500)) +
  scale_fill_gradient(low = "#f2f2f2", high = "#333333",
                      name = NULL,
                      breaks = seq(0, 6000, 1000),
                      guide = NULL) +
  geom_mapa() +
  geom_point(data = square, aes(lon, lat), size = 0.2) +
  geom_point(data = square_d02, aes(lon, lat), size = 0.2, color = "#047994") +
  annotate("text", x = -73.5, y = -40.5, color = "black", label = "d01", size = 3) +
  annotate("text", x = -66, y = -39.5, color = "#047994", label = "d02", size = 3) +
  theme_minimal(base_size = 10) +
  theme(legend.box = "horizontal",
        legend.position = c(0.1, 0.95),
        legend.background = element_rect(fill = "white", color = "white")) 

```


## Resultados

### Sensibilidad a la combinación de canales {#canales}

En este trabajo se generaron 3 experimentos que comparten la misma configuración que los experimentos previamente descriptos pero utilizando distintas combinaciones de canales de ABI para evaluar la sensibilidad a la asimilación de distintas combinaciones de canales. El experimento ABI_ch890_th30 asimila observaciones de los 3 canales, ABI_ch90_th30 asimila observaciones de los canales 9 y 10 (vapor de agua en niveles medios y bajos respectivamente) y ABI_ch0_th30 solo utiliza observaciones del canal 10. Estos 3 experimentos utilizan un thinning de 30 km mientras que un análisis más detallado del thinning se realizará en la siguiente Sección.

(ref:area-canales) Porcentaje de área cubierta por precipitación a 1 hora superior a a) 1 $mmh^{-1}$, b) 5 $mmh^{-1}$ y c) 10 $mmh^{-1}$ a lo largo del tiempo para el campo preliminar de los experimentos SATWND (línea naranja), ABI_ch890_th30 (línea verde), ABI_ch90_th30 (línea violeta) y ABI_ch0_th30 (línea celeste) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre en el dominio de verificación (cuadro celeste en la Figura \@ref(fig:dominio)a). En sombreado y para cada experimento se muestra el área máxima y mínima estimada por el ensamble del análisis. En línea discontinua negra se muestra la estimación de IMERG.

```{r area-canales, fig.cap="(ref:area-canales)", fig.width=6, fig.height=3.5}
imerg <- read_rds(here("data/derived_data/observations/IMERG_1h.rds"))

q <- c(1, 5, 10, 25, 50)
# lon_band = c(-67, -54.5)
# lat_band = c(-38.5, -21)
lon_band = c(-66.5, -61.5)
lat_band = c(-35.5, -29)
imerg <- purrr::map(q, function(f){
  
  imerg %>%
    .[, c("lon", "lat") := wrf_project(x, y, inverse = TRUE, round = FALSE)] %>% 
    .[lon %between% lon_band & lat %between% lat_band] %>%
    .[, .(area = sum(pp_acum > f)/.N), by = .(end_date)] %>%
    .[, umbral := f] %>%
    .[]
}) %>%
  rbindlist()


files <- Sys.glob(here("data/derived_data/analysis_variables/pp_area/E*_ana*_area_acum_*.rds"))

pp_area <- purrr::map(files, readRDS) %>% rbindlist()

pp_area %>% 
  .[umbral %in% c(1, 5, 10) & exp %in% c("E6", "EG2", "EG4", "EG5")] %>% 
  .[, .(mean_area = mean(area),
        max_area = max(area),
        min_area = min(area)), by = .(umbral, date, exp)] %>% 
  ggplot(aes(date, mean_area)) +
  geom_line(data = imerg[umbral %in% c(1, 5, 10)], 
            aes(end_date, area, color = "OBS"),
            linetype = 2) +
  geom_ribbon(aes(ymin = min_area, ymax = max_area, fill = exp), alpha = .1) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(OBS = "grey20", colores_canales),
                     labels = c(E6 = "SATWND",
                                EG2 = "ABI_ch890_th30", EG4 = "ABI_ch90_th30",
                                EG5 = "ABI_ch0_th30", OBS = "OBS")) +
  scale_fill_manual(values = c(OBS = "grey20", colores_canales),
                    guide = NULL,
                    labels = c(E6 = "SATWND",
                               EG2 = "ABI_ch890_th30", EG4 = "ABI_ch90_th30",
                               EG5 = "ABI_ch0_th30", OBS = "OBS")) +
  scale_y_continuous(label = scales::percent_format()) +
  scale_date(20181122000000, 20181123000000, "12 hours", limits = c(ymd_h(2018112200), ymd_h(2018112312))) +
  facet_wrap(~umbral, scales = "free_y", 
             labeller = labeller(umbral = c("1" = "1 mm", "5" = "5 mm", 
                                            "10" = "10 mm", "25" = "25 mm"))) +
  tagger::tag_facets() +
  labs(y = "Área cubierta por precipitación (%)", x = NULL, 
       fill = NULL, color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))

```

La Figura \@ref(fig:area-canales) muestra el porcentaje de área cubierta por precipitación superior a distintos umbrales calculado sobre el pronóstico a 1 hora (campo preliminar) y compara los experimentos para estudiar la sensibilidad al uso de canales con SATWND como experimento control y la estimación de precipitación de IMERG. Al igual que en otros experimentos, los experimentos ABI subestiman la precipitación observada aunque en menor medida que SATWND para el umbral de 1 mm (Figura \@ref(fig:area-canales)a). En algunos casos, logran representar mejor la precipitación que SATWND particularmente para umbrales de precipitación mayores (Figura \@ref(fig:area-canales)b-c) y algunas subregiones como la provincia de Córdoba (no se muestra). Comparando los experimentos ABI entre si, no se observan grandes diferencias y solo en algunos tiempos y/o umbrales ABI_ch890_th30 es mejor que los otros dos. Desde este punto de vista, la asimilación de las observaciones de los canales 8 y 9 en niveles altos y medios no parece aportar mucha información al análisis. Esto podría deberse a que los canales 8 y 9, o sus errores, están muy correlacionados y sería necesario realizar nuevos experimentos que aíslen el impacto de los canales 8 y 9 por separado para determinar si su asimilación independiente no produce mejores resultados que la asimilación del canal 10 o si la asimilación de cualquiera de los 3 canales es suficiente.

Desde el punto de vista del FSS (Figura \@ref(fig:fss-canales)) los experimentos que incluyen observaciones de ABI tienen un impacto positivo mucho mayor que SATWND. Además, los 3 experimentos de ABI tienen un comportamiento similar entre si aunque en este caso ABI_ch890_th30 representa ligeramente mejor la precipitación por encima de 25 mm (Figuras \@ref(fig:fss-canales)b y d). 

Teniendo en cuenta que la asimilación de los 3 canales de ABI en conjunto no degradan el análisis y por el contrario, en algunos casos parecen tener un rendimiento mejor, para los experimentos de este capítulo se asimilaran los 3 canales asociados al vapor del agua al mismo tiempo. 

(ref:fss-canales) FSS calculado sobre la precipitación acumulada a 1 hora en una ventana móvil de 6 horas para umbrales de 1 mm (a y c) y 25 mm (b y d), en escalas de 10 km (a y b) y 100 km (c y d), para el campo preliminar de los experimentos SATWND (línea naranja), ABI_ch890_th30 (línea verde), ABI_ch90_th30 (línea violeta) y ABI_ch0_th30 (línea celeste).

```{r fss-canales, fig.cap="(ref:fss-canales)", fig.width=6, fig.height=5}

files <- c(Sys.glob(here("data/derived_data/FSS/fss_6h_ana_ens_EG[2345]*")),
           Sys.glob(here("data/derived_data/FSS/fss_6h_ana_ens_E[169]*.csv")))

fss <- purrr::map(files, function(f) {
  
  meta <- unglue(basename(f), "fss_6h_ana_ens_{exp}.csv")
  
  fread(f) %>% 
    .[, date := as_datetime(date)] 
}) %>% 
  rbindlist()

fss %>% 
  .[q %in% c(1, 25) & w %in% c(1, 11) & exp %in% c("E6", "EG2", "EG4", "EG5")] %>% 
  .[, ":="(q_label = factor(q, labels = c("1~mm~h^{-1}", "25~mm~h^{-1}")),
           w_label = factor(w, labels = c("10~km", "100~km")))] %>% 
  ggplot(aes(date, fss)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = colores_canales, 
                     labels = c(E6 = "SATWND", 
                                EG2 = "ABI_ch890_th30", EG4 = "ABI_ch90_th30",
                                EG5 = "ABI_ch0_th30")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123120000))) +
  facet_grid(w_label ~ q_label, labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS", linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```


### Sensibilidad al thinning {#thinning}

Existen diversos algoritmos para realizar el thinning, GSI en particular utiliza una metodología simple que consiste en definir una retícula de baja resolución y seleccionar las observaciones en alta resolución de acuerdo la distancia de estás a los puntos de la retícula gruesa y según criterios de calidad descriptos en la Sección \@ref(sat). 

Es importante definir entonces la resolución apropiada para la retícula de baja resolución y estudiar el impacto que esto tiene en el análisis. Para esto se realizaron 2 experimentos ABI_ch890_th30 con un thinning de 30 km, es decir, que lleva las observaciones con una resolución de 2 km a una retícula de 30 km siguiendo trabajos previos [por ejemplo @lee2019; @singh2016] y ABI_ch890_th10 con un thinning de 10 km, igual a la resolución de las simulaciones numéricas. Las 2 resoluciones elegidas tienen en cuenta por un lado el efecto de la correlación entre los errores de las observaciones y la resolución del modelo y por el otro, no perder información importante que podría incluirse en el análisis. En particular, si bien la resolución del modelo en los experimentos es de 10 km, los procesos que puede representar explicitamente son de escala mayor y por lo tanto es posible que asimilar observaciones con una resolución muy cercana al modelo, que contiene información de procesos en escalas no totalmente resueltas por el modelo, no genere el impacto esperado. 

(ref:area-thinning) Porcentaje de área cubierta por precipitación a 1 hora superior a a) 1 $mmh^{-1}$, b) 5 $mmh^{-1}$ y c) 10 $mmh^{-1}$ a lo largo del tiempo para el campo preliminar de los experimentos SATWND (línea naranja), ABI_ch890_th30 (línea verde), y ABI_ch890_th10 (línea celeste) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre en el dominio de verificación (cuadro celeste en la Figura \@ref(fig:dominio)a). En sombreado y para cada experimento se muestra el área máxima y mínima estimada por el ensamble del análisis. En línea continua negra se muestra la estimación de IMERG.

```{r area-thinning, fig.cap="(ref:area-thinning)", fig.width=6, fig.height=3.5}
pp_area %>% 
  .[umbral %in% c(1, 5, 10) & exp %in% c("E6", "EG2", "EG3")] %>% 
  .[, .(mean_area = mean(area),
        max_area = max(area),
        min_area = min(area)), by = .(umbral, date, exp)] %>% 
  ggplot(aes(date, mean_area)) +
  geom_line(data = imerg[umbral %in% c(1, 5, 10)], 
            aes(end_date, area, color = "OBS"),
            linetype = 2) +
  geom_ribbon(aes(ymin = min_area, ymax = max_area, fill = exp), alpha = .1) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(OBS = "grey20", colores_thinning),
                     labels = c(E6 = "SATWND",
                                EG2 = "ABI_ch890_th30", EG3 = "ABI_ch890_th10",
                                OBS = "OBS")) +
  scale_fill_manual(values = c(OBS = "grey20", colores_thinning),
                    guide = NULL,
                    labels = c(E6 = "SATWND",
                               EG2 = "ABI_ch890_th30", EG3 = "ABI_ch890_th10",
                               OBS = "OBS")) +
  scale_y_continuous(label = scales::percent_format()) +
  scale_date(20181122000000, 20181123000000, "12 hours", limits = c(ymd_h(2018112200), ymd_h(2018112312))) +
  facet_wrap(~umbral, scales = "free_y", 
             labeller = labeller(umbral = c("1" = "1 mm", "5" = "5 mm", 
                                            "10" = "10 mm", "25" = "25 mm"))) +
  tagger::tag_facets() +
  labs(y = "Área cubierta por precipitación (%)", x = NULL, 
       fill = NULL, color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))

```

La Figura \@ref(fig:area-thinning) muestra el porcentaje de área cubierta por precipitación superior a distintos umbrales a lo largo del tiempo para los experimentos de thinning y tomando a SATWND como experimento control. En todos los casos ABI_ch890_th10, el experimento con un thinning de 10 km, es el que representa mejor el área cubierta por precipitación, particularmente en las horas de mayor actividad convectiva. Si bien todos los experimentos subestiman la precipitación en comparación con IMERG, ABI_ch890_th10 se acerca a esta estimación sobre todo cuando observamos el valor máximo del ensamble. 

(ref:thinning-fss) FSS calculado sobre la precipitación acumulada a 1 hora en una ventana móvil de 6 horas para umbrales de 1 mm (a y c) y 25 mm (b y d), en escalas de 10 km (a y b) y 100 km (c y d), para el campo preliminar de los experimentos SATWND (línea naranja), ABI_ch890_th30 (línea verde) y ABI_ch890_th10 (línea celeste).

```{r thinning-fss, fig.cap="(ref:thinning-fss)", fig.width=6, fig.height=5}

fss %>% 
  .[q %in% c(1, 25) & w %in% c(1, 11) & exp %in% c("E6", "EG2", "EG3")] %>% 
  .[, ":="(q_label = factor(q, labels = c("1~mm~h^{-1}", "25~mm~h^{-1}")),
           w_label = factor(w, labels = c("10~km", "100~km")))] %>% 
  ggplot(aes(date, fss)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = colores_thinning, 
                     labels = c(E6 = "SATWND", 
                                EG2 = "ABI_ch890_th30", EG3 = "ABI_ch890_th10")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123120000))) +
  facet_grid(w_label ~ q_label, labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS", linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

Desde el punto de vista del FSS y para 1 $mm h^{-1}$ (Figura \@ref(fig:thinning-fss)a y c) ABI_ch890_th10 tiene igual o mejor rendimiento que ABI_ch890_th30 pero ocurre lo contrario para el umbral de 25 $mm h^{-1}$ (Figura \@ref(fig:thinning-fss)b y d). Si bien las diferencias de FSS entre los experimentos de thinning son muy pequeñas, el área cubierta por precipitación y otros análisis como la comparación cualitativa con estimaciones de IMERG y cuantitativa con radiosondeos de RELAMPAGO (no mostrados), resalta la necesidad de utilizar un thinning cercano a la resolución del modelo. Por esta razón la configuración de thinning a utilizar en los experimentos posteriores será de 10 km. 

### Comparación del impacto de asimilar observaciones de ABI con otras fuentes de información

#### Impacto de las observaciones de ABI en el análisis

Al igual que en la Sección \@ref(impacto-analisis), es importante primero entender el impacto de asimilar este conjunto de observaciones en el análisis. En particular, este impacto podría variar en presencia de otros conjuntos de observaciones tales como otras radianzas de satélites polares. Por esta razón se comparará el experimento ABI (previamente ABI_ch890_th10) que asimila los 3 canales de vapor de agua con SATWND donde no se asimilan radianzas y es considerado el experimento control, RAD donde se asimilan radianzas de satélites en órbitas polares y RAD+ABI donde se incluyen todas las observaciones al mismo tiempo. De esta manera podremos evaluar cuál es el impacto de las observaciones de ABI de manera independiente y cuando son asimiladas en conjunto con radianzas de satélites polares.   

En primer lugar calculamos perfiles verticales de temperatura y humedad específica promediando espacialmente la media del ensamble de cada experimento a lo largo de todo el periodo. Posteriormente para evaluar el impacto de la asimilación de las radianzas, calculamos la diferencia entre los perfiles verticales de los experimentos con estas observaciones y SATWND que no las incluye (Figura \@ref(fig:TQ-diff-abi)). Previamente habíamos observado que la asimilación de radianzas de satélites polares producía un calentamiento en niveles bajos (Figura \@ref(fig:TQ-diff-abi)a), la asimilación de observaciones de ABI en este caso produce, además, un enfriamiento en niveles medios y un calentamiento previo al paso del frente frío y desarrollo de la convección durante el 21 de noviembre lo que puede contribuir a incrementar el CAPE y por lo tanto las condiciones para el desarrollo de convección intensa. 

(ref:TQ-diff-abi) Diferencia entre la media del ensamble de los análisis a) y d) RAD-SATWND, b) y e) RAD+ABI-SATWND, y c) y f) ABI-SATWND para los perfiles verticales espacialmente promediados de la temperatura (a, b y c, en $K$) y la humedad específica (d, e y f en $gkg^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis. 

```{r TQ-diff-abi, fig.cap="(ref:TQ-diff-abi)", fig.height=5, fig.width=6, fig.align="left"}
files <- Sys.glob(here("data/derived_data/analysis_variables/perfiles_ana_E*.csv"))

perfiles <- purrr::map(files, function(f) {
  exp <- unglue(basename(f), "perfiles_{run}_{exp}.csv")
  fread(f) %>% 
    .[, exp := exp[[1]][["exp"]]]
  
}) %>% 
  rbindlist() %>% 
  .[exp %in% c("E6", "E9", "E10", "EG3")] %>% 
  .[, date := as_datetime(date)] %>% 
  .[date %between% c(as_datetime("2018-11-20 18:00:00"),
                     as_datetime("2018-11-23 12:00:00"))]


perfiles %>% 
  dcast(lev + date ~ exp, value.var = "T") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)),
                    breaks = c(seq(-1.2, 1, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-1.2, 1, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"), 
                                   labels = scales::number_format()) +
  labs(fill = "Temperatura (K)",
       x = NULL,
       y = "Presión (hPa)") +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6" = "RAD - SATWND",
                                              "E10_E6" = "RAD+ABI - SATWND",
                                              "EG3_E6" = "ABI - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "QVAPOR") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value*1000, fill = stat(level)), 
                    breaks = c(seq(-1.4, 1, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-1.4, 1, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"),
                                   labels = scales::number_format()) +
  
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       y = "Presión (hPa)",
       x = NULL) +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6" = "RAD - SATWND",
                                              "E10_E6" = "RAD+ABI - SATWND",
                                              "EG3_E6" = "ABI - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  plot_layout(ncol = 1, heights = c(1, 1))

```

A las 12 UTC del 22 de noviembre se observa un dipolo en niveles medios y altos que se intensifica con la asimilación de ABI. Esto puede estar relacionado a un menor desarrollo de la convección en ABI a esa hora respecto de SATWND. Pero a partir de las 00 UTC del 23 de noviembre el dipolo se invierte indicando que la convección se intensifica en ABI. Esto también podría explicar el enfriamiento que se observa en superficie debido a una pileta de aire frio más intensa asociada al sistema convectivo.
Al comparar la diferencia entre los perfiles verticales de cada experimento y ERA5 (Figuras \@ref(fig:era5-abi)a-d), vemos que este enfriamiento en niveles medios (Figura \@ref(fig:era5-abi)c-d) produce una mayor diferencia entre los análisis que incluyen ABI y ERA5. 

La diferencia en la humedad específica para los distintos experimentos se muestra en las Figuras \@ref(fig:TQ-diff-abi)d-f, donde vemos que la asimilación de observaciones de ABI produce un humedecimiento por encima de 850 hPa y este efecto es mayor cuando no se asimilan radianzas de satélites polares en el experimento ABI (Figura \@ref(fig:TQ-diff-abi)f). Este humedecimiento junto con los cambios en la temperatura generan un aumento en el CAPE previo al desarrollo de la convección (no se muestra). El secamiento que se observa en los experimentos con radianzas y que ya se observaba previamente en RAD, no es consistente con el desarrollo de la precipitación, que ocurre primero en SATWND consumiendo la humedad presente en la capa límite, por lo que son necesarios nuevos experimentos para entender mejor el mecanismo que genera este efecto. Nuevamente, la diferencia entre los experimentos y ERA5 para esta variable (Figuras \@ref(fig:era5-abi)e-h) muestra que la asimilación de observaciones de ABI aleja los análisis del estado de la atmósfera que representa ERA5.

(ref:era5-abi) Diferencia entre la media del ensamble del análisis de cada experimento y el ERA5 para los perfiles verticales espacialmente promediados de la temperatura del aire (K, a--d), la humedad específica ($g kg^{-1}$, e--h) y el viento meridional ($m s^{-1}$, i--l) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis.

```{r era5-abi, fig.cap="(ref:era5-abi)", fig.width=6, fig.height=8, fig.align="left"}
era <- fread(here("data/derived_data/reanalysis/perfiles_ana_era5.csv"))

era[perfiles, on = c("lev", "date")] %>% 
  .[, exp := factor(exp, levels = c("E6", "E9", "E10", "EG3"))] %>% 
  # .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.T - T, fill = stat(level)), 
                    breaks = seq(-3, 3, 0.5)) +
  geom_contour2(aes(z = i.T - T), color = "white", size = 0.05,
                breaks = seq(-3, 3, 0.5)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E6 = "SATWND", E9 = "RAD", 
                                                              E10 = "RAD+ABI", EG3 = "ABI"))) +
  labs(fill = "Temperatura (K)",
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  era[perfiles, on = c("lev", "date")] %>% 
  .[, exp := factor(exp, levels = c("E6", "E9", "E10", "EG3"))] %>% 
  # .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.QVAPOR - QVAPOR, fill = stat(level)), 
                    breaks = seq(-0.0025, 0.0015, 0.0004)) +
  geom_contour2(aes(z = i.QVAPOR - QVAPOR), color = "white", size = 0.05,
                breaks = seq(-0.0025, 0.0015, 0.0004)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E6 = "SATWND", E9 = "RAD", 
                                                              E10 = "RAD+ABI", EG3 = "ABI"))) +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95), tag_pool = c("e", "f", "g", "h")) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  era[perfiles, on = c("lev", "date")] %>% 
  .[, exp := factor(exp, levels = c("E6", "E9", "E10", "EG3"))] %>% 
  # .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.U - U, fill = stat(level)), 
                    breaks = seq(-4, 5, 0.5)) +
  geom_contour2(aes(z = i.U - U), color = "white", size = 0.05,
                breaks = seq(-4, 5, 0.5)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E6 = "SATWND", E9 = "RAD", 
                                                              E10 = "RAD+ABI", EG3 = "ABI"))) +
  labs(fill = latex2exp::TeX("U ($m$ $s^{-1}$)"),
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95), tag_pool = c("i", "j", "k", "l")) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  era[perfiles, on = c("lev", "date")] %>% 
  .[, exp := factor(exp, levels = c("E6", "E9", "E10", "EG3"))] %>% 
  # .[lev != 1000] %>%
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = i.V - V, fill = stat(level)), 
                    breaks = seq(-4, 7, 1)) +
  geom_contour2(aes(z = i.V - V), color = "white", size = 0.05,
                breaks = seq(-4, 7, 1)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(vars(exp), ncol = 4, labeller = labeller(exp = c(E6 = "SATWND", E9 = "RAD", 
                                                              E10 = "RAD+ABI", EG3 = "ABI"))) +
  labs(fill = latex2exp::TeX("V ($m$ $s^{-1}$)"),
       x = NULL,
       y = "Presión (hPa)") +
  tag_facets(position = list(x = 0.05, y = 0.95), tag_pool = c("m", "n", "ñ", "o")) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        plot.margin = margin(0, 0, 0, 0, "cm"),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  plot_layout(ncol = 1, heights = c(1, 1, 1, 1)) 
```

Las Figuras \@ref(fig:UV-diff-abi) muestran la comparación entre los experimentos que incluyen radianzas y SATWND para las componentes del viento. En este caso las diferencias entre experimentos son más sutiles sin embargo, se aprecia una disminución en el viento zonal, y por lo tanto en la cortante de vertical de viento, por encima de la 700 hPa en los experimentos RAD+ABI y ABI. Esta disminución es más marcada cuando se asimilan radianzas de satélites polares y ABI en conjunto. En este caso la comparación con ERA5 (Figuras \@ref(fig:era5-abi)i-l) muestra que el viento zonal en los experimentos RAD y ABI se acerca más a ERA5 que SATWND que no incluye radianzas o RAD+ABI que asimila todas las radianzas en conjunto. Esto indicaría que los experimentos que incluyen observaciones de ABI representan mejor el flujo zonal para este caso de estudio. En niveles bajos las diferencias son pequeñas y solo se observa un ligero aumento del viento zonal en ABI respecto de SATWND (Figura \@ref(fig:UV-diff-abi)c) asociado al flujo postfrontal. 

Las Figuras \@ref(fig:UV-diff-abi)d-f muestran que los experimentos que asimilan ABI tienen una mayor diferencia respecto de SATWND en el viento meridional en altura a partir del 22 de noviembre a las 12 UTC cuando hay mayor convección en el dominio. Este efecto ya era notorio al analizar el impacto de las radianzas de satélites polares (Figura \@ref(fig:UV-diff)f) que ahora se intensifica cuando se incluyen radianzas de ABI. Estas diferencias están asociadas a un retraso en el desarrollo de la convección en los experimentos que incluyen radianzas respecto de SATWND. Al igual que para el viento meridional, la comparación entre los experimentos y ERA5 (Figuras \@ref(fig:era5-abi)m-o) muestra que, por encima de 500 hPa, los experimentos que incluyen radianzas tienen un viento meridional más parecido a ERA5 que SATWND. Esta mejora es más notoria en ABI donde no se incluyen observaciones de satélites polares.

(ref:UV-diff-abi) Diferencia entre la media del ensamble de los análisis a) y d) RAD-SATWND, b) y e) RAD+ABI-SATWND, y c) y f) ABI-SATWND para los perfiles verticales espacialmente promediados del viento zonal (a, b y c, en $ms^{-1}$) y viento meridional (d, e y f en $ms^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis. Los contornos negros corresponden al viento zonal y meridional para (a,d) RAD, (b,e) RAD+ABI, y (c,f) ABI ya que son los experimentos tienen más observaciones asimiladas en cada panel.

```{r UV-diff-abi, fig.cap="(ref:UV-diff-abi)", fig.height=5, fig.width=6, fig.align="left"}
perfiles %>% 
  dcast(lev + date ~ exp, value.var = "U") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  # geom_contour(aes(z = value)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), breaks = c(seq(-2, 2, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-2, 2, 0.2), Inf)) +
  geom_contour2(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                  melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                breaks = seq(0, 30, 3), size = 0.1, linetype = 1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                    breaks = seq(0, 30, 3), color = "grey30", skip = 1,
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("U ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6"  = "**RAD** - SATWND",
                                              "E10_E6" = "**RAD+ABI** - SATWND",
                                              "EG3_E6" = "**ABI** - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3), 
        strip.text = ggtext::element_markdown()) +
  
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "V") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), 
                    breaks = c(-Inf, seq(-3, 2, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(-Inf, seq(-3, 2, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"),
                                   labels = function(x) JumpBy(x, 2, fill = "")) +
  geom_contour2(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                  melt(id.vars = c("lev", "date")) %>% unique(),
                aes(z = value, linetype = factor(-sign(..level..))),
                breaks = seq(-10, 4, 2), size = 0.1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                    breaks = seq(-10, 4, 2), color = "grey30", skip = 1,
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("V ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_linetype(guide = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6"  = "**RAD** - SATWND",
                                              "E10_E6" = "**RAD+ABI** - SATWND",
                                              "EG3_E6" = "**ABI** - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3),
        strip.text = ggtext::element_markdown()) +
  plot_layout(ncol = 1, heights = c(1, 1))


```



Habiendo analizado el impacto de la asimilación de radianzas en términos generales, la Figura \@ref(fig:sondeos-abi) muestra el RMSE y BIAS calculado al comparar los experimentos con los radiosondeos de RELAMPAGO en la región que se observa la Figura \@ref(fig:dominio)b. Por debajo de 5 km las observaciones de ABI generan una degradación de la temperatura de rocío durante el IOP 7 y de la temperatura en el IOP 8 con aumento tanto del RMSE como del BIAS. Pero se observa una mejora en la representación del viento zonal y meridional en niveles bajos durante el IOP 8. Sin embargo, las mayores diferencias entre SATWND o RAD y los experimentos que incluyen observaciones de ABI se observan por encima de 5 km donde se ubican los máximos de las funciones de peso de los canales de vapor de agua (Figuras \@ref(fig:wf) y \@ref(fig:wf-goes)). La asimilación de observaciones de ABI, en este caso, ayuda a disminuir el RMSE y el BIAS en todas las variables durante el IOP 7 por encima de la capa límite. Este IOP estuvo caracterizado por ser un periodo relativamente estable con cielos mayoritariamente despejados y un flujo del norte en niveles bajos asociado a una advección cálida y húmeda por lo que se asimilaron una mayor cantidad de observaciones de radianzas en comparación con el día siguiente. Esta es una de las razones que podrían explicar el impacto que se observa en todas las variables. Es interesante notar que el error en la humedad durante el IOP 7 mejora tanto en ABI como en RAD+ABI, mientras que RAD es el experimento con mayor RMSE. Esto muestra que la combinación de distintas fuentes de observaciones genera impactos que no son lineales. Durante el IOP 8 las mejoras se observan principalmente en la temperatura y temperatura de punto de rocío y en menor medida en el viento meridional en niveles altos. Finalmente, la característica sobresaliente de este análisis es la mejora sustancial en la representación de la humedad en niveles medios cuando se asimilan observaciones de ABI. Si bien esto es opuesto a lo observado al analizar las diferencias entre ABI y RAD+ABI con ERA5, es importante notar en este punto que la comparación con el reanálisis tiene limitaciones tanto por la resolución del reanálisis como por las observaciones disponibles para su generación.     


(ref:sondeos-abi) RMSE (línea sólida) y BIAS (línea discontinua) de a) la temperatura ($K$), b) la temperatura del punto de rocío ($K$), c) el viento zonal ($m s^{-1}$) y d) el viento meridional ($m s^{-1}$) calculados comparando el análisis de cada experimento con los radiosondeos de RELAMPAGO durante el IOP 7 y el IOP 8. La línea naranja corresponde a SATWND, la línea roja a RAD, RAD+ABI se representa con una línea verde y ABI con una línea turquesa.

```{r sondeos-abi, fig.cap="(ref:sondeos-abi)", fig.width=6, fig.height=4, fig.align="center"}
files <- c(Sys.glob(here("data/derived_data/soundings/sondeo_E[169]*_ana*")),
           Sys.glob(here("data/derived_data/soundings/sondeo_EG3_ana*")))

sondeos <- purrr::map(files, function(f){
  
  out <- fread(f) %>% 
    .[, value := fifelse(variable %in% c("t", "td"), value + 273.15, value)] %>% 
    .[, launch_time := as_datetime(launch_time)]
  
}) %>% 
  rbindlist()

IOP <- tribble(
  ~iop, ~ini, ~end,
  "IOP07", ymd_hms("20181121150000"), ymd_hms("20181121210000"),
  "IOP08", ymd_hms("20181122140000"), ymd_hms("20181122200000")
) %>% setDT()


sondeos %>% 
  .[, iop := fcase(launch_time %between% c(ymd_hms("20181121150000"), ymd_hms("20181121210000")), "IOP07",
                   launch_time %between% c(ymd_hms("20181122140000"), ymd_hms("20181122200000")), "IOP08")] %>% 
  .[variable %in% c("t", "td", "u", "v") & !str_detect(site, "/") & !is.na(fcst_value) & !is.na(iop)] %>% 
  .[site != "Sao Borja, Brazil"] %>% 
  .[, lev := cut_round(alt, c(seq(0, 3000, 500), seq(4000, 21000, 1000)))] %>% 
  .[, ":="(RMSE = mean(sqrt((fcst_value - value)^2), na.rm = TRUE), 
           BIAS = mean(fcst_value - value, na.rm = TRUE)), by = .(iop, exp, variable, lev)] %>% 
  melt(measure.vars = c("RMSE", "BIAS"), variable.name = "estadistico", value.name = "value.est") %>% 
  .[lev <= 15000] %>% 
  .[, ":="(iop_labels = factor(iop, labels = c("IOP~7", "IOP~8")),
           variable_labels = factor(variable, labels = c("Temperatura~(K)", "atop('Temperatura de', 'punto de rocío (K)')", "Viento~U~(m~s^{-1})", "Viento~V~(m~s^{-1})")))] %>% 
  ggplot(aes(lev/1000, value.est)) +
  geom_hline(yintercept = 0, color = "darkgrey") +
  geom_line(aes(color = exp, linetype = estadistico)) +
  scale_color_manual(values = c(colores_exp_abi), labels = c(E6 = "SATWND", E9 = "RAD",
                                                             E10 = "RAD+ABI", EG3 = "ABI")) +
  coord_flip() +
  facet_grid(iop_labels ~ variable_labels, scales = "free_x", 
             labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.95)) +
  labs(y = "BIAS/RMSE", x = "Altura (Km)", color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

##### Impacto en la precipitación

En esta Sección ponemos el foco nuevamente en la precipitación para evaluar la habilidad de los análisis al representar la precipitación acumulada horaria del campo preliminar cuando incluyen o no información de las observaciones de ABI. En primer lugar analizamos el porcentaje de área cubierta por precipitación para distintos umbrales (Figuras \@ref(fig:area-abi)) donde se observa una mejora considerable en la ubicación temporal del área para los experimentos ABI, RAD+ABI y RAD cuando se los compara con SATWND que no asimila radianzas. Sin embargo, la incorporación de observaciones de ABI es lo que genera el mayor impacto en la magnitud del área cubierta por precipitación que aumenta más de un 5% para el umbral de 1 mm en comparación con RAD que solo incluye radianzas de satélites polares. 

En este punto es interesante observar las posibles diferencias entre RAD+ABI y ABI para evaluar el aporte de las radianzas de satélites polares cuando se asimilan en conjunto con radianzas de ABI. Desde el punto de vista del área cubierta por precipitación, las diferencias son pequeñas, RAD+ABI subestima el área previo al máximo de precipitación respecto de ABI y lo contrario ocurre posterior al máximo. Esto podría indicar que las radianzas de ABI son suficientes para representar la precipitación y que la asimilación de radianzas de satélites polares no agrega información, al menos para este caso de estudio. 

(ref:area-abi) Porcentaje de área cubierta por precipitación a 1 hora superior a a) 1 $mmh^{-1}$, b) 5 $mmh^{-1}$ y c) 10 $mmh^{-1}$ a lo largo del tiempo para el campo preliminar de los experimentos SATWND (línea naranja), RAD (línea roja), RAD+ABI (línea verde), y ABI (línea turquesa) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre en el dominio de verificación (cuadro celeste en la Figura \@ref(fig:dominio)a). En sombreado y para cada experimento se muestra el área máxima y mínima estimada por el ensamble. En línea discontinua negra se muestra la estimación de IMERG.

```{r area-abi, fig.cap="(ref:area-abi)", fig.width=6, fig.height=3.5}

pp_area %>% 
  .[umbral %in% c(1, 5, 10) & exp %in% c("E6", "E9", "E10", "EG3")] %>% 
  .[, .(mean_area = mean(area),
        max_area = max(area),
        min_area = min(area)), by = .(umbral, date, exp)] %>% 
  ggplot(aes(date, mean_area)) +
  geom_line(data = imerg[umbral %in% c(1, 5, 10)], 
            aes(end_date, area, color = "OBS"),
            linetype = 2) +
  geom_ribbon(aes(ymin = min_area, ymax = max_area, fill = exp), alpha = .1) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(OBS = "grey20", colores_exp_abi),
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI", OBS = "OBS")) +
  scale_fill_manual(values = c(OBS = "grey20", colores_exp_abi),
                    guide = NULL,
                    labels = c(E6 = "SATWND", E9 = "RAD",
                               E10 = "RAD+ABI", EG3 = "ABI", OBS = "OBS")) +
  scale_y_continuous(label = scales::percent_format()) +
  scale_date(20181122000000, 20181123000000, "12 hours", limits = c(ymd_h(2018112200), ymd_h(2018112312))) +
  facet_wrap(~umbral, scales = "free_y", 
             labeller = labeller(umbral = c("1" = "1 mm", "5" = "5 mm", 
                                            "10" = "10 mm", "25" = "25 mm"))) +
  tagger::tag_facets() +
  labs(y = "Área cubierta por precipitación (%)", x = NULL, 
       fill = NULL, color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))

```

Finalmente, la Figura \@ref(fig:fss-abi) muestra el FSS calculado comparando la estimación de IMERG con cada experimento para distintos umbrales y escalas espaciales, siendo 10 km igual a la resolución del modelo. En todos los casos RAD+ABI y ABI tienen valores mayores de FSS que RAD que ya mostraba mejoras importantes respecto de SATWND para umbrales bajos de precipitación. Solo para el umbral de 25 $mmh^{-1}$ hay un periodo alrededor de las 12 UTC del 22 de noviembre donde todos los análisis que asimilan radianzas no logran representar correctamente la convección y muestran un FSS menor a SATWND. Esto se repite para umbrales de precipitación mayores donde SATWND genera precipitación más intensa pero menos ordenada que ABI durante las primeras 12 horas del 22 de noviembre. Es particularmente notorio si se observan campos de precipitación para cada experimento, donde se ve que el área donde se produce precipitación en ABI es homogénea mientras que en SATWND es discontinua con focos de precipitación aislados (no se muestra). Comparando RAD+ABI y ABI vemos que entre las 06 y las 18 UTC del 22 de noviembre ABI representa la precipitación ligeramente mejor que RAD+ABI, lo que podría indicar que las radianzas de satélites polares tienen un impacto negativo en el análisis. Sin embargo hacia el final del periodo simulado, la relación es opuesta. RAD+ABI tiene valores de FSS mayores que ABI y cercanos a RAD. En la siguiente Sección por lo tanto, se evaluará el impacto de asimilar radianzas en pronósticos a corto plazo y se analizará si el impacto observado en los análisis perdura en los pronósticos y mejora la representación de la convección.

(ref:fss-abi) FSS calculado sobre la precipitación acumulada a 1 hora en una ventana móvil de 6 horas para umbrales de 1 mm (a y c) y 25 mm (b y d), en escalas de 10 km (a y b) y 100 km (c y d), para el campo preliminar de los experimentos SATWND (línea naranja), RAD (línea roja), RAD+ABI (línea verde), y ABI (línea turquesa) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre.

```{r fss-abi, fig.cap="(ref:fss-abi)", fig.width=6, fig.height=5}


fss %>% 
  .[q %in% c(1, 25) & w %in% c(1, 11) & exp %in% c("E6", "E9", "E10", "EG3")] %>% 
  .[, ":="(q_label = factor(q, labels = c("1~mm~h^{-1}", "25~mm~h^{-1}")),
           w_label = factor(w, labels = c("10~km", "100~km")))] %>% 
  ggplot(aes(date, fss)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = colores_exp_abi, 
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123120000))) +
  facet_grid(w_label ~ q_label, labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS", linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```


#### Impacto de las observaciones de ABI en pronósticos determinísticos de alta resolución

A continuación se analizará el impacto de la asimilación de observaciones de ABI en los pronósticos de alta resolución generados a partir de los análisis. También se incluye en la comparación un pronóstico inicializado a partir del análisis global de GFS (generado con el *Global Data Assimilation System*, GDAS), que se utiliza como condición inicial y de borde para la generación de pronósticos operativos en Argentina. 

##### Representación de la temperatura, humedad y viento en los pronósticos

En primer lugar se analiza la representación de la temperatura, humedad y viento en los pronósticos de alta resolución. En este caso se compararan los experimentos RAD, RAD+ABI y ABI con SATWND como experimento control para evaluar el impacto de incorporar las radianzas de satélites polares y de ABI. Posteriormente se comparan los experimentos con observaciones independientes de EMA y radiosondeos de RELAMPAGO.

En primer lugar calculamos perfiles verticales de temperatura y humedad específica promediados espacialmente a partir de los pronósticos determinísticos a 36 hs, con los que generamos la diferencia entre los perfiles verticales de los experimentos con asimilación de radianzas y SATWND. La diferencia en los perfiles promediados de temperatura (Figura \@ref(fig:TQ-diff-fcst)a-c) muestra que las radianzas de satélites polares producen un calentamiento en niveles bajos en comparación con SATWND y esto se repite, aunque en menor medida, con la asimilación de observaciones de ABI. De hecho luego de las 16 UTC del 22 de noviembre, los pronósticos inicializados a partir de RAD+ABI y ABI muestran una capa límite más fría que SATWND lo que podría indicar que desarrollaron un frente frío o una pileta de aire frío asociada a la convección más intensos. En niveles medios y altos las características de las diferencias, si bien de mayor magnitud, reflejan lo observado en los análisis (Figura \@ref(fig:TQ-diff-abi)a-c) por lo que, nuevamente, la información aportada por la asimilación regional de las observaciones produce un impacto prolongado en el pronóstico que determinan el momento de inicio e intensificación de la convección. Respecto de la humedad (Figura \@ref(fig:TQ-diff-fcst)d-f), la comparación con SATWND muestra que la asimilación de observaciones de ABI en las condiciones iniciales luego genera una disminución de la humedad en niveles bajos del pronóstico en comparación con RAD, luego de las 12 UTC del 22 de noviembre y particularmente en las últimas 12 horas de pronóstico. Esta diferencia puede estar asociada al desarrollo de la convección que consume la humedad disponible en capas bajas y que podría ser más intensa en ABI que en RAD.


(ref:TQ-diff-fcst) Diferencia entre los pronósticos determinísticos de a) y d) RAD-SATWND, b) y e) RAD+ABI-SATWND, y c) y f) ABI-SATWND para los perfiles verticales espacialmente promediados de la temperatura (a, b y c, en $K$) y la humedad específica (d, e y f en $g\ kg^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis.

```{r TQ-diff-fcst, fig.cap="(ref:TQ-diff-fcst)", fig.height=5, fig.width=6, fig.align="left"}
perfiles <- Sys.glob(here("data/derived_data/forecast_variables/perfiles_fcst_det*")) %>% 
  map(fread) %>% 
  rbindlist()

perfiles %>% 
  dcast(lev + date ~ exp, value.var = "t") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)),
                    breaks = c(seq(-1.2, 1.2, 0.2), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-1.2, 1.2, 0.2), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3,
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"), 
                                   labels = scales::number_format()) +
  labs(fill = "Temperatura (K)",
       x = NULL,
       y = "Presión (hPa)") +
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6" = "RAD - SATWND",
                                              "E10_E6" = "RAD+ABI - SATWND",
                                              "EG3_E6" = "ABI - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "q") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value*1000, fill = stat(level)), 
                    breaks = c(seq(-3, 1.5, 0.25), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-3, 1.5, 0.25), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"),
                                   labels = scales::number_format()) +
  
  scale_y_level(name = "Presión (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  labs(fill = latex2exp::TeX("Humedad \nespecífica ($g$ $Kg^{-1}$)"),
       y = "Presión (hPa)",
       x = NULL) +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6" = "RAD - SATWND",
                                              "E10_E6" = "RAD+ABI - SATWND",
                                              "EG3_E6" = "ABI - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3)) +
  plot_layout(ncol = 1, heights = c(1, 1))
```

Las diferencias para el viento zonal (Figuras \@ref(fig:UV-diff-fcst)a-c) muestran que la asimilación de observaciones de satélites polares y ABI en conjunto produce el mayor impacto en los pronósticos respecto de SATWND en todos los niveles pero en particular por encima de la capa límite. En cuanto al viento meridional (Figuras \@ref(fig:UV-diff-fcst)d-f), se observa que los pronósticos con condiciones iniciales que incluyen observaciones de ABI son más cercanos a SATWND que previamente había demostrado ser el pronóstico con mejor representación de la precipitación (ver Sección \@ref(prono-impacto)). Por otro lado, la evolución del viento meridional a lo largo del tiempo para cada pronóstico (contornos en la Figuras \@ref(fig:UV-diff-fcst)d-f) muestra que los experimentos ABI y RAD+ABI tiene un flujo en altura por detrás del sistema convectivo es más intenso (Figuras \@ref(fig:UV-diff-fcst)e-f) que RAD. El dipolo que se observa en la diferencia entre los pronósticos, nuevamente indica que la convección que genera RAD o ABI está retrasada en el tiempo respecto de SATWND. Sin embargo, las diferencias entre ABI-SATWND disminuyen respecto de RAD-SATWND, lo que podría ser el resultado del uso de condiciones iniciales con mayor información del entorno convectivo.

Si bien todo lo anterior corresponde a una comparación relativa al pronóstico SATWND, en líneas generales el signo del impacto observado en los pronósticos inicializados a partir de los análisis que asimilan radianzas (Figuras \@ref(fig:TQ-diff-fcst) y \@ref(fig:UV-diff-fcst)) coincide con lo observado en la evaluación del impacto de los análisis (Figuras \@ref(fig:TQ-diff-abi) y \@ref(fig:UV-diff-abi)).

(ref:UV-diff-fcst) Diferencia entre los pronósticos determinísticos de a) y d) RAD-SATWND, b) y e) RAD+ABI-SATWND, y c) y f) ABI-SATWND para los perfiles verticales espacialmente promediados del viento zonal (a, b y c, en $ms^{-1}$) y viento meridional (d, e y f en $ms^{-1}$) calculados sobre el dominio interior (recuadro rojo en la Figura \@ref(fig:dominio)a) para cada ciclo de análisis. Los contornos negros corresponden al viento zonal y meridional para (a,d) RAD, (b,e) RAD+ABI, y (c,f) ABI.

```{r UV-diff-fcst, fig.cap="(ref:UV-diff-fcst)", fig.height=5, fig.width=6, fig.align="left"}
perfiles %>% 
  dcast(lev + date ~ exp, value.var = "u") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  # geom_contour(aes(z = value)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), breaks = c(seq(-4, 4, 0.5), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black")) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(seq(-4, 4, 0.5), Inf)) +
  geom_contour2(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                  melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                breaks = seq(0, 30, 3), size = 0.1, linetype = 1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                    breaks = seq(0, 30, 3), color = "grey30", skip = 1,
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("U ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6"  = "**RAD** - SATWND",
                                              "E10_E6" = "**RAD+ABI** - SATWND",
                                              "EG3_E6" = "**ABI** - SATWND"))) +
  tag_facets(position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3), 
        strip.text = ggtext::element_markdown()) +
  
  
  perfiles %>% 
  dcast(lev + date ~ exp, value.var = "v") %>% 
  .[, ":="(E9_E6 = E9 - E6,
           E10_E6 = E10 - E6,
           EG3_E6 = EG3 - E6)] %>% 
  melt(measure.vars = c("E9_E6", "E10_E6", "EG3_E6")) %>% 
  ggplot(aes(date, lev)) +
  geom_contour_fill(aes(z = value, fill = stat(level)), 
                    breaks = c(-Inf, seq(-6, 6, 0.5), Inf)) +
  geom_contour2(aes(z = value), color = "white", size = 0.05,
                breaks = c(-Inf, seq(-6, 6, 0.5), Inf)) +
  scale_fill_divergent_discretised(guide = guide_colorsteps(barwidth = 25,
                                                            barheight = 0.3, 
                                                            title.position = "left", 
                                                            title.vjust = 1,
                                                            frame.colour = "black"),
                                   labels = function(x) JumpBy(x, 2, fill = "")) +
  geom_contour2(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                  melt(id.vars = c("lev", "date")) %>% unique(),
                aes(z = value, linetype = factor(-sign(..level..))),
                breaks = seq(-24, 10, 2), size = 0.1, color = "grey30") +
  geom_text_contour(data = ~.x[, .(lev, date, E9_E6 = E9, E10_E6 = E10, EG3_E6 = EG3)] %>%
                      melt(id.vars = c("lev", "date")) %>% unique(), aes(z = value),
                    breaks = seq(-24, 10, 4), color = "grey30", skip = 1,
                    size = 2, stroke = 0.1) +
  labs(fill = latex2exp::TeX("V ($m$ $s^{-1}$)"),
       x = NULL,
       y = NULL) +
  scale_linetype(guide = NULL) +
  scale_y_level(name = "Pressure (hPa)", breaks = c(1000, 850, 750, 500, 300, 200, 100)) +
  scale_date(ini = 20181121000000, end = 20181123110000, break_bin = "12 hours") +
  facet_wrap(~variable, ncol = 3,
             labeller = labeller(variable = c("E9_E6"  = "**RAD** - SATWND",
                                              "E10_E6" = "**RAD+ABI** - SATWND",
                                              "EG3_E6" = "**ABI** - SATWND"))) +
  tag_facets(tag_pool = c("d", "e", "f"), position = list(x = 0.05, y = 0.96)) +
  theme_minimal(base_size = 8) +
  theme(legend.position = "bottom",
        tagger.panel.tag.text = element_text(size = 8),
        panel.ontop = TRUE,
        panel.grid = element_line(linetype = 3),
        strip.text = ggtext::element_markdown()) +
  plot_layout(ncol = 1, heights = c(1, 1))


```

En la comparación de los pronósticos con observaciones independientes se incluye también el pronóstico determinístico inicializado a partir del análisis de GFS. En la Figura \@ref(fig:sondeos-abi-fcst) se muestran los perfiles de RMSE y BIAS calculados comparando los pronósticos con los radiosondeos de RELAMPAGO. La Figura \@ref(fig:sondeos-abi-fcst)a muestra que los errores asociados a la temperatura en niveles bajos son mayores en los experimentos que incluyen radianzas de satélites polares en comparación con GFS o ABI. En niveles medios y altos GFS presenta el mayor RMSE y BIAS en la temperatura, lo que pone en evidencia el impacto positivo de la asimilación de observaciones a escala regional. Esto mismo ocurre para la temperatura de punto de rocío (Figura \@ref(fig:sondeos-abi-fcst)b) que muestra que el RMSE y BIAS de GFS es hasta 5$^\circ$K mayor que para el resto de los pronósticos en niveles medios. El RMSE del viento zonal también muestra que GFS es el pronóstico con mayor error alrededor de 5 km mientras que por encima de los 10 km el error en RAD aumenta por sobre los otros (Figura \@ref(fig:sondeos-abi-fcst)c). Para el viento meridional, el RMSE de todos los pronósticos tiene un comportamiento similar aunque en algunos niveles (alrededor de 5 y 11 km) SATWND muestra un aumento en el error. En lineas generales ninguno de los pronósticos inicializados a partir de análisis con asimilación regional es consistentemente mejor para todas las variables analizadas. 

(ref:sondeos-abi-fcst) RMSE (línea sólida) y BIAS (línea discontinua) de a) la temperatura ($K$), b) la temperatura del punto de rocío ($K$), c) el viento zonal ($m\ s^{-1}$) y d) el viento meridional ($m\ s^{-1}$) calculados comparando los pronósticos determinísticos de cada experimento con los radiosondeos de RELAMPAGO. La línea naranja corresponde a SATWND, la línea roja a RAD, RAD+ABI se representa con una línea verde, ABI con una línea turquesa y en negro se muestra la comparación con el pronóstico inicializado con GFS.

```{r sondeos-abi-fcst, fig.cap="(ref:sondeos-abi-fcst)", fig.width=6, fig.height=4, fig.align="center"}
files <- Sys.glob(here("data/derived_data/soundings/20181122000000_det/sondeo_*"))

sondeos <- purrr::map(files, function(f){
  
  fread(f) %>% 
    .[, value := fifelse(variable %in% c("t", "td"), value + 273.15, value)] %>% 
    .[, launch_time := as_datetime(launch_time)]
  
}) %>% 
  rbindlist()

IOP <- tribble(
  ~iop, ~ini, ~end,
  "IOP07", ymd_hms("20181121150000"), ymd_hms("20181121210000"),
  "IOP08", ymd_hms("20181122140000"), ymd_hms("20181122200000")
) %>% setDT()


sondeos %>% 
  .[variable %in% c("t", "td", "u", "v") & !str_detect(site, "/") & !is.na(fcst_value)] %>% 
  .[site != "Sao Borja, Brazil"] %>% 
  .[, lev := cut_round(alt, c(seq(0, 3000, 500), seq(4000, 21000, 1000)))] %>% 
  .[, ":="(RMSE = mean(sqrt((fcst_value - value)^2), na.rm = TRUE), 
           BIAS = mean(fcst_value - value, na.rm = TRUE)), by = .(exp, variable, lev)] %>% 
  melt(measure.vars = c("RMSE", "BIAS"), variable.name = "estadistico", value.name = "value.est") %>% 
  .[lev <= 15000] %>% 
  .[, ":="(variable_labels = factor(variable, labels = c("Temperatura~(K)", "atop('Temperatura de', 'punto de rocío (K)')", "Viento~U~(m~s^{-1})", "Viento~V~(m~s^{-1})")))] %>% 
  ggplot(aes(lev/1000, value.est)) +
  geom_hline(yintercept = 0, color = "darkgrey") +
  geom_line(aes(color = exp, linetype = estadistico)) +
  scale_color_manual(values = c(colores_exp_abi, GFS = "grey20"),  
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI")) +
  coord_flip() +
  facet_grid(. ~ variable_labels, scales = "free_x", 
             labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.95)) +
  labs(y = "BIAS/RMSE", x = "Altura (Km)", color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

El RMSE y BIAS calculado comparando los pronósticos con las observaciones de EMA se muestra en la Figura \@ref(fig:aut-rmse-fcst) para todo el periodo pronosticado. Si bien en general todos los pronósticos tienen errores similares, se observa que los pronósticos inicializados a partir de análisis con asimilación regional tienen errores grandes al inicio del pronóstico que luego disminuyen rápidamente en la primera hora. Esto posiblemente esté asociado a que el sistema de asimilación solo actualiza las variables 3D del modelo durante la asimilación y no así la temperatura o humedad a 2 metros que se utiliza en esta verificación. Otro factor que podría influir es el uso de la media del ensamble como condición inicial, tema que se discutirá más adelante. Por otro lado el BIAS de SATWND es consistentemente menor en las primeras 10 horas de pronóstico posiblemente gracias a la asimilación de observaciones de EMA que mejoran la representación de la humedad en las condiciones iniciales. Para la temperatura, GFS mantiene un error menor al resto de los experimentos durante las últimas 12 horas de pronóstico aunque el comportamiento general de los errores no difiere del resto. Cómo fue analizado previamente, ninguno de los pronósticos inicializados a partir de análisis con asimilación regional es consistentemente mejor para todas las variables analizadas hasta ahora. 

(ref:aut-rmse-fcst) Evolución del RMSE y BIAS de a) humedad específica ($g kg{^-1}$), b) temperatura ($K$), c) el viento zonal ($m\ s^{-1}$) y d) el viento meridional ($m s^{-1}$) calculados comparando los pronósticos determinísticos de cada experimento con las observaciones de EMA. La línea naranja corresponde a SATWND, la línea roja a RAD, RAD+ABI se representa con una línea verde, ABI con una línea turquesa y en negro se muestra el pronóstico inicializado con GFS.

```{r aut-rmse-fcst, fig.cap="(ref:aut-rmse-fcst)", fig.align="center", fig.width=6, fig.height=3}

files_obs <- Sys.glob(here("data/derived_data/forecast_variables/interp_conv/interp_conv_*"))

ema <- map(files_obs, read_rds) %>% 
  rbindlist()

ema %>% 
  .[usage.flag == 1] %>% 
  melt(id.vars = c("lon", "lat", "var", "date", "obs"), measure.vars = c("E10", "E6", "E9", "EG3", "GFS"),
       variable.name = "exp") %>% 
  .[, value := fifelse(var == "q", value*1000, value)] %>%
  .[, obs := fifelse(var == "q", obs*1000, obs)] %>%
  .[, .(RMSE = sqrt(mean((value - obs)^2, na.rm = TRUE)),
        BIAS = mean(value - obs, na.rm = TRUE)), by = .(var, date, exp)] %>%
  melt(measure.vars = c("RMSE", "BIAS")) %>% 
  .[, ":="(estadistico_label = factor(variable),
           var_label = factor(var, labels = c("atop('Humedad', 'específica (g Kg'^-1*')')",
                                              "Temperatura~(K)",
                                              "Viento~U~(m~s^{-1})",
                                              "Viento~V~(m~s^{-1})")))] %>% 
  ggplot(aes(date, value)) +
  geom_hline(data = hline, aes(yintercept = value), color = "darkgrey") +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(colores_exp_abi, GFS = "grey20"),  
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI")) +
  # facet_grid(variable~var, scales = "free") +
  scale_date(20181122000000, 20181123000000, "12 hours", 
             limits = c(ymd_h(2018112200), ymd_h(2018112312))) +
  labs(x = NULL, y = NULL, linetype = NULL, color = NULL) +
  facet_grid(estadistico_label ~ var_label, scales = "free_y",
             labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        aspect.ratio = 1/1.5,
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```


##### Representación de la convección y precipitación

Para analizar la representación de la precipitación se compararon los pronósticos determinísticos con las estimaciones de IMERG para calcular el FSS y el porcentaje de área cubierta por precipitación. En cuanto a la convección, se utilizó el modelo de transferencia radiativa CRTM para calcular la temperatura de brillo asociada al canal 13 (infrarrojo) de ABI a partir de los pronósticos y así compararlos cualitativamente con las observaciones de GOES-16. Además se calculó el FSS para distintos umbrales de temperatura de brillo para hacer una comparación cuantitativa del desarrollo de la convección observada por el satélite GOES-16.

El porcentaje de área cubierta por precipitación (Figura \@ref(fig:area-abi-fcst)) muestra que los pronósticos subestiman la ocurrencia de precipitación aún cuando la resolución del modelo es alta y representa la convección de manera explicita. Si bien los pronósticos inicializados a partir de GFS y SATWND tienen un mayor porcentaje de área cubierta por precipitación mayor a 1 mm y en algunos tiempos también para precipitación mayor a 5 mm, los pronósticos que incluyen información de las observaciones de ABI representan muy bien el máximo observado para todos los umbrales. También es interesante notar que SATWND tiene un comportamiento similar al resto de los pronósticos y apenas se observan indicios de que la precipitación se adelante en el tiempo. Esto podría deberse a que los pronósticos inicializados con observaciones de ABI contribuyen a la iniciación temprana de la precipitación disminuyendo las diferencias entre SATWND y ABI. GFS es el pronóstico que genera precipitación más rápidamente, ya que el área cubierta en todos los umbrales es mayor durante las primeras 6 a 12 horas. Sin embargo, hacia el final del periodo pronosticado GFS es el que muestra un porcentaje de área cubierta menor en comparación con el resto de los pronósticos.

(ref:area-abi-fcst) Porcentaje de área cubierta por precipitación a 1 hora superior a a) 1 $mmh^{-1}$, b) 5 $mmh^{-1}$ y c) 10 $mmh^{-1}$ a lo largo del tiempo para el campo preliminar de los experimentos SATWND (línea naranja), RAD (línea roja), RAD+ABI (línea verde), ABI (línea turquesa) y GFS (línea negra continua) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre en el dominio *d02* (Figura \@ref(fig:dominio-det)). En sombreado y para cada experimento se muestra el área máxima y mínima estimada por el ensamble del análisis. En línea discontinua negra se muestra la estimación de IMERG.

```{r area-abi-fcst, fig.cap="(ref:area-abi-fcst)", fig.width=6, fig.height=3.5}
imerg <- read_rds(here("data/derived_data/observations/IMERG_1h.rds"))

q <- c(1, 5, 10, 25, 50)

lon_band = c(-67.7, -53.7)
lat_band = c(-40.5, -21.3)

imerg <- purrr::map(q, function(f){
  
  imerg %>%
    .[, c("lon", "lat") := wrf_project(x, y, inverse = TRUE, round = FALSE)] %>% 
    .[lon %between% lon_band & lat %between% lat_band] %>%
    .[, .(area = sum(pp_acum > f)/.N), by = .(end_date)] %>%
    .[, umbral := f] %>%
    .[]
}) %>%
  rbindlist()

files <- Sys.glob(here("data/derived_data/analysis_variables/pp_area/*_fcst_det*area*"))

pp_area <- purrr::map(files, function(f) {
  
  meta <- unglue(f, "{exp}_fcst_det_{ini}_area_acum_1h.rds")
  
  readRDS(f) %>% 
    .[, ini_date := meta[[1]][["ini"]]]
  
}) %>% rbindlist()

pp_area %>% 
  .[umbral %in% c(1, 5, 10) & exp %in% c("E6", "E9", "E10", "EG3", "GFS")] %>% 
  .[ini_date == "2018112200"] %>% 
  .[, .(mean_area = mean(area),
        max_area = max(area),
        min_area = min(area)), by = .(umbral, date, exp)] %>% 
  ggplot(aes(date, mean_area)) +
  geom_line(data = imerg[umbral %in% c(1, 5, 10)], 
            aes(end_date, area, color = "OBS"),
            linetype = 2) +
  # geom_ribbon(aes(ymin = min_area, ymax = max_area, fill = exp), alpha = .1) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(OBS = "grey20", colores_exp_abi, GFS = "grey20"),
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI",
                                GFS = "GFS", OBS = "OBS")) +
  scale_y_continuous(label = scales::percent_format()) +
  scale_date(20181122000000, 20181123000000, "12 hours", limits = c(ymd_h(2018112200), ymd_h(2018112312))) +
  facet_wrap(~umbral, scales = "free_y", 
             labeller = labeller(umbral = c("1" = "1 mm", "5" = "5 mm", 
                                            "10" = "10 mm", "25" = "25 mm"))) +
  tagger::tag_facets() +
  labs(y = "Área cubierta por precipitación (%)", x = NULL, 
       fill = NULL, color = NULL, linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))

```

El FSS calculado comparando los pronósticos de precipitación con la estimación de IMERG (Figura \@ref(fig:fss-abi-fcst)) muestra nuevamente que GFS representa de manera más precisa la precipitación al comienzo del pronóstico, seguido por SATWND. Sin embargo con el paso de las horas de pronóstico, los experimentos con asimilación de radianzas pasan a tener valores de FSS mayores que GFS y esto se mantiene hasta el final del periodo pronosticado. Esto indica, para este caso de estudio, que si bien los pronósticos inicializados a partir de análisis con asimilación regional demoran en producir precipitación, luego mantienen la generación de la misma por más tiempo. Para esta métrica no se ven diferencias apreciables entre los distintos umbrales o escalas espaciales. 

(ref:fss-abi-fcst) FSS calculado sobre la precipitación acumulada en una ventana móvil de 6 horas para umbrales de 1 mm (a y d), 5 mm (b y e) y 25 mm (c y f), en escalas de 10 km (a-c) y 100 km (d-f), para los pronósticos determinísticos SATWND (línea naranja), RAD (línea roja), RAD+ABI (línea verde), ABI (línea turquesa) y GFS (línea negra) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre.

```{r fss-abi-fcst, fig.cap="(ref:fss-abi-fcst)", fig.width=5, fig.height=4}
files <- Sys.glob(here("data/derived_data/FSS/fss_6h_fcst_det_2018112200*d02.csv"))


fss_det <- map(files, function(f){
  
  meta <- unglue(basename(f), "fss_6h_fcst_det_{ini_date}_{exp}_d02.csv")
  
  fread(f) %>% 
    .[, ini_date := meta[[1]][["ini_date"]]]
}) %>% 
  rbindlist()

fss_det[q %in% c(1, 5, 10) & w %in% c(1, 5) & exp %in% c("E6", "E9", "E10", "EG3", "GFS")] %>% 
  .[, ":="(q_label = factor(q, labels = c("1~mm~h^{-1}",  "5~mm~h^{-1}", "25~mm~h^{-1}")),
           w_label = factor(w, labels = c("2~km", "10~km")))] %>% 
  ggplot(aes(date, fss)) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(colores_exp_abi, GFS = "grey20"),  
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123100000))) +
  facet_grid(w_label ~ q_label, labeller = label_parsed) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS", linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))
```

Por otro lado, la Figura \@ref(fig:tb) muestra la temperatura de brillo simulada para el canal 13 de GOES-16 en el tiempo inicial del pronóstico a las 00 UTC del 22 de noviembre y luego a las 12 UTC del mismo día cuando el sistema convectivo comienza a desarrollarse y a las 00 UTC del 23 de noviembre cuando el sistema convectivo ya se encuentra en el norte de Argentina y presenta mayor intensidad. A las 00 UTC del 22 de noviembre los pronósticos, particularmente SATWND, desarrollan nubosidad estratiforme y de poco espesor en gran parte del dominio que no está presente en las observaciones (Figura \@ref(fig:tb)a-f). Sin embargo los pronósticos inicializados a partir de RAD+ABI y ABI ubica mejor la nubosidad y la generación de cúmulos con desarrollo vertical al norte de la provincia de Salta. Además al incluir información sobre las regiones con cielo despejado, gracias a la asimilación de radianzas de ABI, también representa mejor las zonas despejadas al este del dominio. En este sentido, el pronóstico inicializado a partir de ABI logra reducir en algunos tiempos, la convección espuria y a veces también la precipitación que genera SATWND (no se muestra). Debido a que el GFS no se inicializa con información sobre el contenido de agua de nube no muestra nubosidad para el primer tiempo de pronóstico (Figura \@ref(fig:tb)e).

A las 12 UTC, cuando la convección se encuentra cruzando el dominio, todos los pronósticos ubican la nubosidad en la región correcta aunque menos organizada y asociadas a temperaturas de tope más altas que lo observado. De todos los pronósticos SATWND es el que muestra valores de temperatura de brillo más bajos (Figura \@ref(fig:tb)g). Por otro lado, la nubosidad que se observa sobre Chaco también está presente en los pronósticos inicializados a partir del análisis con asimilación observaciones de ABI (Figuras \@ref(fig:tb)i y j). Sin embargo no logran generar la nubosidad sobre Uruguay que si se muestra en el pronóstico de GFS (Figuras \@ref(fig:tb)k). 

Finalmente a las 00 UTC del 23 de noviembre la convección se extiende sobre todo el noreste del dominio con temperaturas de brillo por debajo de los -85$^\circ$C asociada a los topes nubosos más altos (Figura \@ref(fig:tb)r). Si bien todos los pronósticos logran representar las características principales del sistema convectivo la temperatura de brillo es ligeramente mayor y solo en algunas regiones llega a valores por debajo de -85$^\circ$C. En este caso GFS es el único pronóstico que representa con mayor detalle la región con menores temperaturas de brillo dividida en 2 subregiones como se observa en la imagen de brillo de GOES-16.

(ref:tb) Temperatura de brillo simulada por CRTM a partir de los pronósticos determinísticos a las 00 y 12 UTC del 22 de noviembre y las 00 UTC del 23 de noviembre. Las observaciones de temperatura de brillo del canal 13 de GOES-16 se muestran en la ultima columna (f, l, r). 

```{r tb, fig.cap="(ref:tb)", fig.width=6.5, fig.height=5, fig.align="left"}
exps <- c("E6", "E9", "E10", "EG3", "GFS")
dates <- c(rep("20181122000000", 5), rep("20181122120000", 5), rep("20181123000000", 5))

files_goes <- Sys.glob(here(paste0("data/derived_data/observations/", exps, "/G16_", dates, ".nc")))

goes_exp <- map(files_goes, function(f) {
  
  meta <- unglue(f, "/home/paola.corrales/tesis_doctorado/data/derived_data/observations/{exp}/G16_{date}.nc")
  ReadNetCDF(f, vars = c("CH13", lon = "XLONG", lat = "XLAT")) %>% 
    .[, ":="(exp = meta[[1]][["exp"]],
             date = ymd_hms(meta[[1]][["date"]]))]
  
}) %>% rbindlist()

goes <- map(here(paste0("data/derived_data/observations/", "interp_goes_", c("20181122000000", "20181122120000", "20181123000000"), ".rds")), read_rds) %>% 
  rbindlist() 


goes <- goes_exp[, .(south_north, west_east, lon, lat)] %>% 
  unique() %>% 
  .[goes, on = .NATURAL] %>% 
  .[, exp_factor := factor("OBS", levels = c("E6", "E9", "E10", "EG3", "GFS", "OBS"))]

goes_exp %>% 
  .[, date := factor(date)] %>% 
  .[, exp_factor := factor(exp, levels = c("E6", "E9", "E10", "EG3", "GFS", "OBS"))] %>% 
  .[date != "2018-11-23 12:00:00"] %>% 
  ggplot(aes(lon, lat)) +
  scattermore::geom_scattermore(aes(color = CH13 - 273.15)) +
  scattermore::geom_scattermore(data = goes, aes(color = tb_goes)) +
  scale_color_topes(guide = guide_colorbar(barwidth = 15,
                                           barheight = 0.5,
                                           order = 99,
                                           frame.colour = "black")) +
  geom_mapa() +
  coord_sf(ylim = c(-41, -21), xlim = c(-68.5, -53)) +
  facet_grid(date ~ exp_factor, labeller = labeller(exp_factor = c("E6" = "SATWND", "E9" = "RAD", 
                                                                   "E10" = "RAD+ABI", "EG3" = "ABI", 
                                                                   "GFS" = "GFS", "OBS" = "OBS"),
                                                    date = c("2018-11-22 00:00:00" = "00 UTC, 22 de Nov",
                                                             "2018-11-22 12:00:00" = "12 UTC, 22 de Nov",
                                                             "2018-11-23 00:00:00" = "00 UTC, 23 de Nov"))) +
  tagger::tag_facets(position = list(x = 0.04, y = 0.95)) +
  labs(color = latex2exp::TeX("TB ($^{o}C$)")) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom")
```

Para analizar la temperatura de brillo desde un punto de vista cuantitativo, la Figura \@ref(fig:fss-tb) muestra el FSS calculado para umbrales de temperatura de brillo de -30$^\circ$C, asociado al borde externo de la nubosidad, -70$^\circ$C que corresponde a la región de colores rojizos y -85$^\circ$C asociado a los topes nubosos. Al igual que con la precipitación, GFS tiene una mejor representación de la temperatura de brillo desde el comienzo del periodo pronosticado pero luego es igual o peor que los pronósticos inicilizados a partir de los análisis. Para el umbral de -30$^\circ$C no se ven diferencias apreciables entre los experimentos y tampoco en las distintas escalas. Sin embargo es destacable que los valores de FSS son cercanos a 1. Al observar el umbral de -70$^\circ$C, vemos que solo en escalas de 200 km hay valores de FSS  para las temperaturas de brillo más bajas lo que indica que las regiones donde la nubosidad es más alta y por lo tanto más fría en los pronósticos no coincide exactamente con las regiones observadas durante las primeras horas de pronóstico. Finalmente, la temperatura de brillo por debajo de -85$^\circ$C, es la más difícil de pronosticar porque solo se observa en regiones muy pequeñas dentro de la nube convectiva donde se encuentran las corrientes ascendentes más intensas. Sumado a esto, cabe mencionar las dificultades del modelo de transferencia radiativa para simular estas regiones de baja temperatura. En este caso los pronósticos simulan temperaturas de brillo de esta magnitud pero en regiones muy limitadas y que no coinciden exactamente con las observaciones. Solo al mirar la escala de 200 km se ve que tanto ABI como GFS logran representar topes nubosos intensos. Es importante notar que si bien el FSS de los pronósticos inicializados a partir de la asimilación regional es alto, luego cae rápidamente debido a una disminución en el contenido de condensados en el modelo. Esto ocurre porque si bien los pronósticos son inicializados usando la información de condensados de la media del ensamble del análisis, que surgen de promediar la información de todos los miembros, estos no son campos físicos válidos. El modelo los descarta durante el ininio de la similación y genera los condensados desde cero, requiriendo tiempo de spin up. Este es solo un problema para lo pronósticos determinísicos generados a partir de la media del ensamble y deberían aplicarse técnicas para inicializar estas variables y pruebas adicionales para evitar el spin up del modelo.

(ref:fss-tb) FSS calculado sobre la temperatura de brillo simulada para umbrales de -30$^\circ$C (a, d, g), -70$^\circ$C (b, e, h) y -85$^\circ$C (c, f, i), en escalas de 2 km (a-c), 20 km (d-f) y 200 km (g-i), para los pronósticos determinísticos SATWND (línea naranja), RAD (línea roja), RAD+ABI (línea verde), ABI (línea turquesa) y GFS (línea negra) entre las 00 UTC del 22 de noviembre y las 12 UTC del 23 de noviembre.

```{r fss-tb, fig.cap="(ref:fss-tb)", fig.width=5, fig.height=4}

files <- Sys.glob(here("data/derived_data/FSS/fss_fcst_det_tb_*"))

fss <- map(files, fread) %>% rbindlist()

fss %>% 
  .[q %in% c(30, 70, 85) & w %in% c(1, 11, 101)] %>% 
  ggplot(aes(date, fss)) +
  geom_line(aes(color = exp)) +
  scale_color_manual(values = c(colores_exp_abi, GFS = "grey20"),  
                     labels = c(E6 = "SATWND", E9 = "RAD",
                                E10 = "RAD+ABI", EG3 = "ABI")) +
  scale_date(20181122000000, 20181123000000, "12 hours") +
  coord_cartesian(xlim = c(ymd_hms(20181122000000), ymd_hms(20181123100000))) +
  # facet_grid(w_label ~ q_label, labeller = label_parsed) +
  facet_grid(w ~ q, labeller = labeller(w = c("1" = "2 km",
                                              "11" = "20 km",
                                              "101" = "200 km"),
                                        q = c("30" = "-30 ºC",
                                              "70" = "-70 ºC",
                                              "85" = "-85 ºC"))) +
  tagger::tag_facets(position = list(x = 0.05, y = 0.90)) +
  labs(color = NULL, x = NULL, y = "FSS", linetype = NULL) +
  theme_minimal(base_size = 9) +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "#fbfbfb", color = NA))



```



## Conclusiones

En este capítulo evaluamos el impacto que produce la asimilación de observaciones de los canales de vapor de agua del sensor ABI a bordo del satélite geoestacionario GOES-16 tanto en el análisis como en pronósticos determinísticos de alta resolución a corto plazo. La asimilación de estas observaciones es de particular interés en Sudamérica por su alta resolución espacial y frecuencia temporal en comparación con otras fuentes de observación como los sensores a bordo de satélites polares. Además, al momento de la escritura de esta tesis, no se conocen trabajos que utilicen observaciones de sensor ABI en la región, por lo que estos resultados podrían contribuir al desarrollo de los sistemas de pronósticos en Sudamérica.

Para asimilar estas observaciones fue primero necesario incorporar las rutinas necesarias en el sistema de asimilación GSI. Por esto y por que se trata de observaciones nuevas, se generaron experimentos de sensibilidad a la combinación de canales de vapor de agua y al thinning y se estudió la distribución de los errores de las observaciones para determinar la mejor configuración del sistema de asimilación. En este trabajo, entonces, se utilizaron los 3 canales de vapor de agua que aportan información en niveles altos (~350 hPa), medios (~400 hPa) y bajos (~600 hPa) de la atmósfera, con una resolución de 10 km y sin aplicar ningún tipo de corrección de bias. 

En primer lugar, de la comparación de los análisis con y sin asimilación de radianzas podemos ver que estas produce impactos en todos los niveles de la atmósfera. En algunos casos estos impactos son directos, es decir, correcciones introducidas a las variables del modelo durante la asimilación, y en otros casos indirectos, es decir causados por procesos físicos que se desencadenan a partir de las correcciones que introducen las observaciones. Si bien los experimentos incrementan ligeramente el bias cálido y seco presente en niveles bajos del modelo, también generan cambios en la temperatura y humedad en niveles medios que son consistentes con un aumento de la inestabilidad previo al desarrollo del sistema convectivo. Los perfiles de las componentes de viento muestran que la iniciación de la convección en los experimentos que asimilan ABI tiene un retraso con respecto del experimento control sin asimilación de radianzas. Sin embargo, la comparación con estimaciones de precipitación muestran que la asimilación de observaciones de ABI mejora sustancialmente el área cubierta por precipitación y la representación general de la precipitación particularmente para umbrales bajos. Esto último podría deberse en parte a una mejor representación de la humedad en niveles medios y altos, algo que resulta evidente al comparar los distintos experimentos con los radiosondeos de RELAMPAGO. 

Los pronósticos determinísticos en alta resolución inicializados a partir de la media del análisis de cada experimento y que resuelven los procesos convectivos explicitamente nos permiten analizar el desarrollo de la convección en más detalle. Es notorio que las condiciones iniciales tienen un impacto a lo largo de todo el pronóstico e influyen en el momento de la iniciación de la convección y su intensidad. Si bien analizando el área cubierta por precipitación o la comparación con las estimaciones de precipitación, todos los pronósticos tiene un comportamiento similar, al compararlos con un pronóstico sin asimilación regional, vemos mejoras tanto en la representación de la temperatura y humedad en niveles medios y altos como en la precipitación en las últimas horas de pronóstico. Esto es de alguna manera sorprendente, ya que es esperable que el origen de las mejoras esté asociado a las condiciones iniciales de los pronósticos. Si bien no hay una razón establecida para este comportamiento, es posible que problemas en la inicialización de los condensados en los pronósticos regionales tenga un impacto negativo en las primeras horas de simulación.

Un análisis cualitativo del desarrollo de la convección a partir de la temperatura de brillo del canal de 10.3 $\mu m$ de GOES-16 muestra que los pronósticos con condiciones iniciales que tienen asimilación de ABI representan muy bien el desarrollo convectivo y también  las regiones del dominio donde hay cielos despejados, eliminando en muchos casos celdas convectivas generadas por el modelo en otros experimentos. Por otro lado, del análisis cuantitativo de la temperatura de brillo surge que solo el pronóstico con observaciones de ABI que no incluye radianzas de satélites polares y el pronóstico sin asimilación regional, logran representar los valores de más bajos de temperatura asociados a las ascendentes intensas en los sistemas convectivos. 


Finalmente, de la comparación entre los experimentos es importante evaluar el impacto de las observaciones de ABI en presencia o no de otras radianzas de satélites polares. Para este caso de estudio se observó que, en general, la asimilación de radianzas de satélites polares no produjo un impacto por encima de la asimilación de radianzas de ABI cuando se las asimilaba en conjunto. Además, en algunos casos los resultados asociados al experimento que incluye solo radianzas de ABI fueron mejores que en los experimentos donde se incluían también radianzas de satélites polares. Es posible entonces, que la asimilación de radianzas de satélites polares, con todos los desafíos que implica no sea necesaria si se asimilan radianzas de ABI para generar una buena representación de procesos convectivos, aunque esta hipótesis necesita confirmarse mediante el análisis de un numero mayor de casos.