01-chap1.Rmd

<!--
This is for including Chapter 1.  Notice that it's also good practice to name your chunk.  This will help you debug potential issues as you knit.  The chunk above is called intro and the one below is called chapter1.  Feel free to change the name of the Rmd file as you wish, but don't forget to change it here from chap1.Rmd.
-->

<!--
The {#rmd-data} text after the chapter declaration will allow us to link throughout the document back to the beginning of Chapter 1.  These labels will automatically be generated (if not specified) by changing the spaces to hyphens and capital letters to lowercase.  Look for the reference to this label at the beginning of Chapter 2.
-->


```{r load_pkgs1, echo=FALSE, message=FALSE, warning=FALSE}
# List of packages required for this analysis
  pkg <- c("dplyr", "sf", "ggplot2","rnaturalearth", "rnaturalearthdata", "ggspatial", "kableExtra", "ncdf4", "stars", "magick", "RcmdrMisc", 
           "knitr", "ggrepel", "grid", "gridExtra", "cowplot", "xts", "bookdown", "lubridate", "devtools")
# Check if packages are not installed and assign the
# names of the packages not installed to the variable new.pkg
  new.pkg <- pkg[!(pkg %in% installed.packages())]
# If there are any packages in the list that aren't installed,
# install them
  if (length(new.pkg))
    install.packages(new.pkg, repos = "http://cran.rstudio.com")
# Load packages (thesisdown will load all of the packages as well)
  library(thesisdown)
  library(dplyr)
  library(sf)
  library(ggplot2)
  library(rnaturalearth)
  library(rnaturalearthdata)
  library(ggspatial)
  library(knitr)
  library(kableExtra)
  library(ncdf4)
  library(stars)
  library(magick)
  library(RcmdrMisc)
  library(ggrepel)
  library(grid)
  library(gridExtra)
  library(cowplot)
  library(xts)
  library(lubridate)
```


```{r loadstations, include=FALSE}
#Load IDEAM and ISD Stations
con1 = src_postgres(dbname = "winddata", host = "localhost", port = 5432, user = "user1", password = "user1")

#Get Ideam Stations Table
originalfields4 = c("objectid", "codigo1", "nombre", "latitud", "longitud", "categoria")
originalfields4 = paste(originalfields4, collapse= ", ", sep = "")
query4 = paste("select", originalfields4, "from ideam_all_stations", "where inpqrs2 = 'YES'", sep=" ")
ideam_stations = as_tibble(tbl(con1, sql(query4)))
Encoding(ideam_stations$categoria) <- "UTF-8"
Encoding(ideam_stations$nombre) <- "UTF-8"

originalfields3 = c("id", "usaf", "station_name", "latitud", "longitud")
originalfields3 = paste(originalfields3, collapse= ", ", sep = "")
query3 = paste("select", originalfields3, "from isd_all_stations where usaf_isd_dataua != ''", sep=" ")
isd_stations = as_tibble(tbl(con1, sql(query3)))
```

# Data {#rmd-data}

Input data was obtained from three different sources (a) Institute of Hydrology, Meteorology and Environmental Studies of Colombia IDEAM <http://www.ideam.gov.co>, (b) Integrated Surface Database ISD <https://www.ncdc.noaa.gov/isd>, and (c) climate reanalysis ERA5 <https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5>. Information about datasets, variables, and units is found in tables \@ref(tab:tabledatasources1), \@ref(tab:tabledatasources2), and \@ref(tab:tabledatasources3) respectively.

```{r tabledatasources1, echo= F, echo=FALSE, message=FALSE, warning=FALSE, results="asis"}
data_sources_thesis_summary1 <-read.csv("./data/data_sources_thesis_summary1.csv",header=TRUE, stringsAsFactors = F)

kable(data_sources_thesis_summary1, align=rep('l', 3),
      #col.names = c("Name[Code]", "Latitude", "Longitude"),
      caption = "Datasets description",
      caption.short = "Datasets Description",
      longtable = TRUE,
      booktabs = TRUE) %>%
      row_spec(0, align = "l") %>%
      column_spec(2,width = "1.2in") %>%
      column_spec(3,width = "3.8in") %>%
kable_styling(font_size = 10, latex_options="scale_down")

```

```{r tabledatasources2, echo= F, echo=FALSE, message=FALSE, warning=FALSE, results="asis"}
data_sources_thesis_summary2 <-read.csv("./data/data_sources_thesis_summary2.csv",header=TRUE, stringsAsFactors = F)

kable(data_sources_thesis_summary2, align=rep('l', 3),
      #col.names = c("Name[Code]", "Latitude", "Longitude"),
      caption = "Variables in all datasets",
      caption.short = "Datasets Variables",
      longtable = TRUE,
      booktabs = TRUE) %>%
      row_spec(0, align = "l") %>%
      column_spec(2,width = "1.2in") %>%
      column_spec(3,width = "3.5in") %>%
kable_styling(font_size = 10, latex_options="scale_down")

```


```{r tabledatasources3, echo= F, echo=FALSE, message=FALSE, warning=FALSE, results="asis"}
data_sources_thesis_summary3 <-read.csv("./data/data_sources_thesis_summary3.csv",header=TRUE, stringsAsFactors = F)

kable(data_sources_thesis_summary3, align=rep('l', 4),
      #col.names = c("Name[Code]", "Latitude", "Longitude"),
      caption = "Variables Units and Time",
      caption.short = "Variables Units and Time",
      longtable = TRUE,
      booktabs = TRUE) %>%
      row_spec(0, align = "l") %>%
      column_spec(2,width = "0.8in") %>%
      column_spec(3,width = "3in") %>%
kable_styling(font_size = 10, latex_options="scale_down")
```

Ideal data source to create extreme wind speeds maps should be field observed data from IDEAM, but the IDEAM data have the following deficiencies: 

1. There are not enough number of stations around the study area to represent all the local wind variability in a huge country with multiple variety of climates and changing thermal floors.

2. There are uncertainties related to the way anemometers are registering data, then comparison with other data sources are needed to be able to do appropriate data standardization, needed as a prerequisite for the analysis.

3. There is no time continuity in the registration of IDEAM data. Historical time series are different and variable in each station.

4. The different wind related variables provided by IDEAM are not well documented and their comparative values for identical time periods are not consistent with each other.

Other logistical difficulties in accessing IDEAM data source are also highlighted:

1. It was not possible to obtain the wind hourly mean variable.

2. It was impossible for them to calculate and deliver a representative gust velocity, for example _5-minutes gust_.

3. There was also no access to complete raw data (instantaneous measurements every minute that they claim to have), to be able to calculate from them the variables required for the study.

Importance of ISD database for this study is based on the fact that post-processed ISD database has wind extreme values, and it was used to create extreme wind maps for United States. ISD allows comparison with IDEAM records to take better decisions in order to conduct data standardization. Despite that ERA5 data are not observed data (but forecast), its main advantage is their resolution (0.25 square decimal degrees) and availability.

## IDEAM

Historical observed wind speeds from `r length(ideam_stations$codigo1)` stations in Colombia are managed by the official environmental authority IDEAM. Table \@ref(tab:tableideamstations) shows a sample of ten IDEAM stations. Figure \@ref(fig:plotideamstations) shows a map of IDEAM stations. Figure \@ref(fig:plotoneideamstation) shows data for IDEAM station "21205791".

```{r tableideamstations, echo=FALSE, message=FALSE, warning=FALSE, results="asis"}
tableideam =  ideam_stations[c(81,120,200,57,179, 22,35,150,170,185), c("nombre", "latitud", "longitud")]
tableideam$latitud = round(tableideam$latitud, 2)
tableideam$longitud = round(tableideam$longitud, 2)

#tableideamstations = 
kable(tableideam, 
      col.names = c("Name[Code]", "Latitude", "Longitude"),
      caption = "IDEAM Stations Sample",
      caption.short = "IDEAM Stations Sample",
      longtable = TRUE,
      booktabs = TRUE)  %>%
kable_styling(font_size = 10, latex_options="scale_down")

```

```{r sfideam, echo=FALSE, message=FALSE, warning=FALSE}
#Create simple features from Ideam Stations
ideam_stations = st_as_sf(ideam_stations, coords = c("longitud", "latitud"), crs = 4326)
```

```{r plotideamstations, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="IDEAM Stations. Colombia"}
#Plot IDEAM Stations
theme_set(theme_bw())
world <- ne_countries(scale = "medium", returnclass = "sf")
world_points<- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))


colombia = world_points$name == "Colombia"
panama = world_points$name == "Panama"
peru= world_points$name == "Peru"
brazil= world_points$name == "Brazil"
venezuela= world_points$name == "Venezuela"
ecuador= world_points$name == "Ecuador"

ggplot(data = world) + 
  geom_sf(fill= "antiquewhite") + 
  geom_text(data= world_points[venezuela,],aes(x=-67, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-79.2, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-78.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-67, y=-2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 13, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = ideam_stations, size=1, aes(shape=categoria, color=categoria), show.legend = "point") + 
  scale_color_discrete(name = 'Category', labels = c("Agrometeorological", "Ordinary Climatic", "Main Climatic", "Mareographic", 
   "Special Meteorological", "Main Synoptic")) +
  scale_shape_discrete(name = 'Category', labels = c("Agrometeorological", "Ordinary Climatic", "Main Climatic", "Mareographic", 
   "Special Meteorological", "Main Synoptic")) +  
  annotation_scale(location = "bl", width_hint = 0.5) + 
  annotation_north_arrow(location = "br", which_north = "true", pad_x = unit(0.05, "in"), pad_y = unit(0.05, "in"), 
   style = north_arrow_fancy_orienteering) + 
  coord_sf(xlim = c(-82.1, -63.8), ylim = c(-7.5, 15.5), expand = FALSE) +
  xlab("Longitude") + 
  ylab("Latitude") + 
  ggtitle("IDEAM Stations") + 
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"))
```

```{r plotoneideamstation, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Time Series of IDEAM Station ELDORADO CATAM - AUT", fig.width=4.5, fig.height=3.5}
#Plote time series - one IDEAM Station
originalfields = c("21205791")
newfields = paste ("X", originalfields, sep="")
originalfields = paste('"', originalfields, '"', sep = "")
newfields = paste('"', newfields, '"', sep = "")

fiedls_query = paste(originalfields, "as", newfields, sep = " ")
fiedls_query = c(paste('"', "mydatetime", '"', sep = ""), fiedls_query)
fiedls_query = paste (fiedls_query, "", sep= "", collapse=", ")

wherestring = c("21205791")
wherestring = paste('"', wherestring, '"', sep = "")
wherestring = paste(wherestring, "IS NOT NULL", sep = " ")
wherestring = paste(wherestring, collapse = " OR " , sep = " ")
query = paste("select", fiedls_query, "from ideam_vvmx_60", "where", wherestring, sep=" ")

all_vvmx_aut_60 = as_tibble(tbl(con1, sql(query)))
timestamp_all_vvmx_aut_60 <- as.POSIXct(as_tibble(select(all_vvmx_aut_60, mydatetime))$mydatetime,format="%Y-%m-%d %H:%M:%S", tz="UTC")

statideam_xts = na.omit(xts(x=select(all_vvmx_aut_60, "X21205791"), order.by = timestamp_all_vvmx_aut_60))

plot.xts(statideam_xts, main = "Station ID: 21205791\nWind Velocity [m/s]", major.ticks="years", format.labels = "%b-%d\n%Y", legend.loc = "top", 
  col="green", cex.main=0.3, cex=0.4, cex.axis=0.9, mar = c(2.5,1,0,1), oma = c(0,0,0,0))
```

```{r plotoneideamstationacf, include=FALSE, fig.cap="IDEAM Station ACF"}
# par(oma = c(0,0,0,0))
# myacf = acf(statideam_xts[,1][!is.na(statideam_xts[,1])], plot=FALSE)
# titleacf = paste("Autocorrelation Function - ACF", paste("IDEAM", "21205791", sep = ":"), sep = "\n")
# plot(myacf, main = titleacf)
```

```{r plotoneideamstationpacf, include=FALSE, fig.cap="IDEAM Station PACF"}
# par(oma = c(0,0,0,0))
# mypacf = pacf(statideam_xts[,1][!is.na(statideam_xts[,1])], plot=FALSE)
# titlepacf = paste("Partial Autocorrelation Function - PACF", paste("IDEAM", "21205791", sep = ":"), sep = "\n")
# plot(mypacf, main = titlepacf)
```

## ISD

ISD is a database with environmental variables, among them extreme wind speeds. ISD has data for the whole planet, and is based on observed data at meteorological stations in each country, which means that for Colombia is based on IDEAM data. Main advantage is data availability at neighbor countries and specialized post-processing made by NOAA's National Centers for Environmental Information NCEI in United States, which facilitates its use. Table \@ref(tab:tableisdstations) shows a sample of ten ISD stations. Figures \@ref(fig:plotisdstations) and \@ref(fig:plotoneisdstation) shows a map of ISD stations and data from ISD station "802590".


```{r tableisdstations, echo=FALSE, message=FALSE, warning=FALSE, results="asis"}
tableisd =  isd_stations[c(81,20,29,57,60, 30,50,55,58,3), c("usaf", "station_name", "latitud", "longitud")]
tableisd$latitud = round(tableisd$latitud, 2)
tableisd$longitud = round(tableisd$longitud, 2)

#tableideamstations = 
kable(tableisd, 
      col.names = c("Code", "Name", "Latitude", "Longitude"),
      caption = "ISD Stations Sample",
      caption.short = "ISD Stations Sample",
      longtable = TRUE,
      booktabs = TRUE) %>%
kable_styling(font_size = 10, latex_options="scale_down")
```

```{r isdsf, echo=FALSE, message=FALSE, warning=FALSE}
#Create simple features from ISD stations
isd_stations = st_as_sf(isd_stations, coords = c("longitud", "latitud"), crs = 4326)
```


```{r plotisdstations, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="ISD Stations. Colombia and Surroundings countries"}
#Plot ISD Stations
theme_set(theme_bw())
world <- ne_countries(scale = "medium", returnclass = "sf")
world_points<- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))

colombia = world_points$name == "Colombia"
panama = world_points$name == "Panama"
peru= world_points$name == "Peru"
brazil= world_points$name == "Brazil"
venezuela= world_points$name == "Venezuela"
ecuador= world_points$name == "Ecuador"

ggplot(data = world) + 
  geom_sf(fill= "antiquewhite") + 
  geom_text(data= world_points[venezuela,],aes(x=-67, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-79.2, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-78.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-67, y=-2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 13, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = isd_stations, size=1, aes(color= "ISD Stations"), shape=2, show.legend = "point") + 
  scale_color_manual(values = c("ISD Stations" = "black"), name="") +
  annotation_scale(location = "bl", width_hint = 0.5) + 
  annotation_north_arrow(location = "br", which_north = "true", pad_x = unit(0.05, "in"), pad_y = unit(0.05, "in"), 
   style = north_arrow_fancy_orienteering) + 
  coord_sf(xlim = c(-82.1, -63.8), ylim = c(-7.5, 15.5), expand = FALSE) +
  xlab("Longitude") + 
  ylab("Latitude") + 
  ggtitle("Integrated Surface Database - ISD") + 
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"))
```

```{r plotoneisdstation, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="ISD Station ALFONSO BONILLA ARAGON INTL - Time Series", fig.width=4.5, fig.height=3.3}
#Plot - ISD Station
originalfields1 = c("802590")
newfields1 = paste ("X", originalfields1, sep="")
originalfields1 = paste('"', originalfields1, '"', sep = "")
newfields1 = paste('"', newfields1, '"', sep = "")
fiedls_query1 = paste(originalfields1, "as", newfields1, sep = " ")
fiedls_query1 = c(paste('"', "mydatetime", '"', sep = ""), fiedls_query1)
fiedls_query1 = paste (fiedls_query1, "", sep= "", collapse=", ")

wherestring1 = c("802590")
wherestring1 = paste('"', wherestring1, '"', sep = "")
wherestring1 = paste(wherestring1, "IS NOT NULL", sep = " ")
wherestring1 = paste(wherestring1, collapse = " OR " , sep = " ")
query1 = paste("select", fiedls_query1, "from isd_lite_unstack", "where", wherestring1, sep=" ")

isdlite = as_tibble(tbl(con1, sql(query1)))

timestamp_isdlite <- as.POSIXct(as_tibble(select(isdlite, mydatetime))$mydatetime,format="%Y-%m-%d %H:%M:%S", tz="UTC")

statisd_xts = na.omit(xts(x=select(isdlite, "X802590"), order.by = timestamp_isdlite))

plot.xts(statisd_xts, main = "Station ID: 802590\nWind Velocity [m/s]", major.ticks="years", format.labels = "%b-%d\n%Y", legend.loc = "top", 
 col="green", cex.main=0.3, cex=0.4, cex.axis=0.9, mar = c(2.5,1,0,1), oma = c(0,0,0,0))
```

```{r plotoneisdstationacf, include=FALSE, fig.cap="ISD Station ACF"}
# par(oma = c(0,0,0,0))
# myacf = acf(statisd_xts[,1][!is.na(statisd_xts[,1])], plot=FALSE)
# titleacf = paste("Autocorrelation Function - ACF", paste("ISD", "802590", sep = ":"), sep = "\n")
# plot(myacf, main = titleacf)
```

```{r plotoneisdstationpacf, include=FALSE, fig.cap="IDEAM Station PACF"}
# par(oma = c(0,0,0,0))
# mypacf = pacf(statisd_xts[,1][!is.na(statisd_xts[,1])], plot=FALSE)
# titlepacf = paste("Partial Autocorrelation Function - PACF", paste("ISD", "802590", sep = ":"), sep = "\n")
# plot(mypacf, main = titlepacf)
```


## ERA5


```{r loadera5, include=FALSE}
#Load ERA5 NetCDF dataset - variable fg10
 ncname <- "outfile_nc4c_zip9"
 filename <- paste("./data/", ncname, ".nc", sep = "")
# ncin <- nc_open(filename)
# lon <- ncvar_get(ncin, "longitude")
# nlon = dim(lon)
# lat <- ncvar_get(ncin, "latitude")
# nlat = dim(lat)
# ntime <-  dim(ncvar_get(ncin, "time"))
# variablename <- "fg10"
# fg10.units <- ncatt_get(ncin, variablename, "units")
# fg10.units
# lonlat.unstack <- expand.grid(lon=as.numeric(lon), lat=as.numeric(lat))

#Load using stars
lon <- read_ncdf(filename, var = c("longitude"))
nlon = dim(lon)
lat <- read_ncdf(filename, var = c("latitude"))
nlat = dim(lat)
ntime <-  dim(read_ncdf(filename, var = c("time")))

variablename <- "fg10"
fg10_1 = read_ncdf(filename, var = c(variablename), ncsub = cbind(start = c(1,1,1), count = c(1,1,1)))
fg10.units = units(fg10_1$fg10)

lonlat.unstack <- expand.grid(lon=as.numeric(lon$longitude), lat=as.numeric(lat$latitude))
```


```{r era5sf, include=FALSE}
#Create ERA5 centers (point with lat, lon, and value, as cell index)
era5colpoints = st_as_sf(lonlat.unstack, coords=1:2, crs=st_crs(4326))
era5colpoints$value = 1:(nlon*nlat)
```

ERA5 is forecast reanalysis data processed by the _European Centre for Medium-Range Weather Forecasts_ ECMWF with wind speeds time series in square cells of 0.25 decimal degrees covering the whole planet. It was extracted a raster of `r (nlat)` rows by `r nlon` columns in format NetCDF. Cell centers represent ERA5 stations, with IDs from 1 (lon=-79, lat=12.5) to 3381 (lon=-67, lat=-4.5). Map in figure \@ref(fig:plotera5stations) shows ERA5 stations.

```{r era5bbox, include=FALSE}
#Define stars object to match with ERA5 bounding box.
#Cell centers of stars object, need to be same cell centers of ERA5
pointsbbox = st_bbox(era5colpoints)
cellsize = lonlat.unstack$lon[2]- lonlat.unstack$lon[1]
mybbox = st_bbox(c(pointsbbox$xmin - (cellsize/2), pointsbbox$xmax + (cellsize/2), pointsbbox$ymax + (cellsize/2), pointsbbox$ymin - (cellsize/2)), 
 crs = st_crs(4326))
era5colraster.st = st_rasterize(era5colpoints, st_as_stars(mybbox, nx = nlon, ny = nlat, values = era5colpoints$value))
```


```{r era5pols, message=FALSE, warning=FALSE, include=FALSE}
#Load ERA5 polygon vectors, representing cells of ERA5
file_era5_sf_pol = "./data/era5grid_left_right_pol.shp"
era5_4326_sf_pol = st_read(dsn=file_era5_sf_pol)
pts <- do.call(rbind, st_centroid(st_geometry(era5_4326_sf_pol)))
x = pts[,1]
y = pts[,2]
era5_4326_sf_pol$x = x
era5_4326_sf_pol$y = y
era5_4326_sf_pol_filter_corners = era5_4326_sf_pol %>% filter(DN %in% c(1, 49, 3333, 3381))
era5_4326_sf_pol_filter_corners_left = era5_4326_sf_pol %>% filter(DN %in% c(1, 3333))
era5_4326_sf_pol_filter_corners_right = era5_4326_sf_pol %>% filter(DN %in% c(49, 3381))
```


```{r plotera5stations, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="ERA5 Cells and Stations (Cells Centers)", fig.height=4.05}
#Plot ERA5 Stations
theme_set(theme_bw())
world <- ne_countries(scale = "medium", returnclass = "sf")
world_points<- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))

colombia = world_points$name == "Colombia"
panama = world_points$name == "Panama"
peru= world_points$name == "Peru"
brazil= world_points$name == "Brazil"
venezuela= world_points$name == "Venezuela"
ecuador= world_points$name == "Ecuador"

big = ggplot(data = world) + 
  geom_sf(fill= "antiquewhite") +
  geom_sf(data=era5colpoints, size=0.1, aes(color = "Stations"), shape=".",  show.legend = "point")+
  scale_color_manual(values = c("Stations" = "black"), name="ERA5", guide = guide_legend(override.aes = list(fill= c(NA), linetype = c("blank"), 
   shape = c(".")))) +
  geom_sf(data = era5_4326_sf_pol_filter_corners, color = "black", aes(fill="Cells"), size=0.1, alpha=1, show.legend = "polygon") +
  scale_fill_manual(values = c("Cells" = NA), name="", guide = guide_legend(override.aes = list(fill = c(NA), shape = c(NA), size=0.1))) +
  geom_rect(mapping=aes(xmin=-79.252968100, xmax=-78.247031900, ymin=11.832846362, ymax=12.667153638), color="red", alpha=0, size=0.1) +
  geom_rect(mapping=aes(xmin=-67.752968100, xmax=-66.747031900, ymin=11.832846362, ymax=12.667153638), color="red", alpha=0, size=0.1) +
  geom_rect(mapping=aes(xmin=-79.258632089, xmax=-78.241367911, ymin=-4.671851259, ymax=-3.828148741), color="red", alpha=0, size=0.1) +
  geom_rect(mapping=aes(xmin=-67.758632089, xmax=-66.741367911, ymin=-4.671851259, ymax=-3.828148741), color="red", alpha=0, size=0.1) +
  geom_text(data= world_points[venezuela,],aes(x=-66.3, y=8.5, label=name), color = "darkblue", fontface = "bold", size=2, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-79.7, y=9.2, label=name), color = "darkblue", fontface = "bold", size=2, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-79.5, y=-1.5, label=name), color = "darkblue", fontface = "bold", size=2, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75.5, y=-5.2, label=name), color = "darkblue", fontface = "bold", size=2, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-66, y=-2, label=name), color = "darkblue", fontface = "bold", size=2, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 14, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 2) + 
  annotate(geom = "text", x = -80, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 2) +
  geom_text_repel(data = era5_4326_sf_pol_filter_corners_left, size=2, aes(x=x, y=y, label = DN), direction="y", segment.size=0.1, 
   segment.color= "grey50", color="grey50", nudge_x=-1, hjust=1, box.padding=0.1) +
  geom_text_repel(data = era5_4326_sf_pol_filter_corners_right, size=2, aes(x=x, y=y, label = DN), direction="y", segment.size=0.1, 
   segment.color= "grey50", color="grey50", nudge_x=+1, hjust=0, box.padding=0.1) +
  coord_sf(xlim = c(-81.1, -64.8), ylim = c(-5, 13), expand = FALSE) +
  xlab("") + 
  ylab("") + 
  ggtitle("ERA5 Reanalysis - Forecast") +
  theme(plot.title = element_text(size=8)) +
  theme(axis.text.x= element_text(size=7)) + 
  theme(axis.text.y= element_text(size=7)) +
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.1)) +
  theme(panel.background = element_rect(fill = "aliceblue")) +
  theme(legend.title = element_text(size=8)) +
  theme(legend.text=element_text(size=8)) +
  theme(legend.key.size = unit(0.5,"line")) +
  theme(plot.margin=unit(c(0,0,0,0),"cm")) +
  theme(axis.text.x = element_text(margin =  margin(t =2, b = -10))) + 
  theme(axis.text.y = element_text(margin =  margin(r =2, l = -10)))
  
corner1lt = ggplot(data = world) + 
  geom_sf(fill= "antiquewhite",  size=0.1) + 
  geom_sf(data = era5_4326_sf_pol, colour="black", fill=NA, size=0.1) +
  geom_text(data= world_points[venezuela,],aes(x=-66.5, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-80.5, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-79.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-68, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[colombia,],aes(x=-71, y=4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 14, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80.5, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = st_cast(world, "MULTILINESTRING"),  size=0.1)+
  geom_sf_text(data = era5_4326_sf_pol, aes(label = DN), size=2) +
  coord_sf(xlim = c(-79.252968100, -78.247031900), ylim = c(11.832846362, 12.667153638), expand = FALSE) +
  xlab("") + 
  ylab("") + 
  ggtitle("") +
  theme(panel.grid = element_blank()) + 
  theme(panel.background = element_rect(fill = "aliceblue")) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
  theme(axis.ticks = element_blank()) +
  theme(plot.margin=grid::unit(c(0,0.2,0,0),"cm")) + 
  theme(panel.border = element_rect(colour = "red"))+
  theme(axis.ticks.length=unit(0, "null")) +
  theme(axis.title.x=element_blank()) +
  theme(axis.title.y=element_blank()) +
  theme(plot.title = element_blank())
  
corner2rt = ggplot(data = world) + 
  geom_sf(fill= "antiquewhite",  size=0.1) + 
  geom_sf(data = era5_4326_sf_pol, colour="black", fill=NA, size=0.1) + 
  geom_text(data= world_points[venezuela,],aes(x=-66.5, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-80.5, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-79.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-68, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[colombia,],aes(x=-71, y=4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 14, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80.5, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = st_cast(world, "MULTILINESTRING"),  size=0.1)+
  geom_sf_text(data = era5_4326_sf_pol, aes(label = DN), size=2) +
  coord_sf(xlim = c(-67.752968100, -66.747031900), ylim = c(11.832846362, 12.667153638), expand = FALSE) +
  xlab("") + 
  ylab("") + 
  ggtitle("") +
  theme(panel.grid = element_blank()) +
  theme(panel.background = element_rect(fill = "aliceblue")) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
  theme(axis.ticks = element_blank()) +
  theme(plot.margin=grid::unit(c(0,0,0,0.2),"cm")) + 
  theme(panel.border = element_rect(colour = "red"))+
  theme(axis.ticks.length=unit(0, "null")) +
  theme(axis.title.x=element_blank()) +
  theme(axis.title.y=element_blank()) +
  theme(plot.title = element_blank())

corner3lb = ggplot(data = world) + 
  geom_sf(fill= "antiquewhite",  size=0.1) + 
  geom_sf(data = era5_4326_sf_pol, colour="black", fill=NA, size=0.1) + 
  geom_text(data= world_points[venezuela,],aes(x=-66.5, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-80.5, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-79.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-68, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[colombia,],aes(x=-71, y=4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 14, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80.5, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = st_cast(world, "MULTILINESTRING"),  size=0.1)+
  geom_sf_text(data = era5_4326_sf_pol, aes(label = DN), size=2) +
  coord_sf(xlim = c(-79.258632089, -78.241367911), ylim = c(-4.671851259, -3.828148741), expand = FALSE) +
  xlab("") + 
  ylab("") + 
  ggtitle("") +
  theme(panel.background = element_rect(fill = "aliceblue")) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
  theme(axis.ticks = element_blank()) +
  theme(plot.margin=grid::unit(c(0,0.2,0,0),"cm")) + 
  theme(panel.border = element_rect(colour = "red"))+
  theme(axis.ticks.length=unit(0, "null")) +
  theme(axis.title.x=element_blank()) +
  theme(axis.title.y=element_blank()) +
  theme(plot.title = element_blank())

corner4rb = ggplot(data = world) + 
  geom_sf(fill= "antiquewhite",  size=0.1) + 
  geom_sf(data = era5_4326_sf_pol, colour="black", fill=NA, size=0.1) + 
  geom_text(data= world_points[venezuela,],aes(x=-66.5, y=8.5, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[panama,],aes(x=-80.5, y=9.2, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) + 
  geom_text(data= world_points[ecuador,],aes(x=-79.2, y=-1, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[peru,],aes(x=-75, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[brazil,],aes(x=-68, y=-6, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  geom_text(data= world_points[colombia,],aes(x=-71, y=4, label=name), color = "darkblue", fontface = "bold", size=3, check_overlap = FALSE) +
  annotate(geom = "text", x = -77.5, y = 14, label = "Caribbean\nSea", fontface = "italic", color = "grey22", size = 4) + 
  annotate(geom = "text", x = -80.5, y = 5, label = "Pacific\nSea", fontface = "italic", color = "grey22", size = 4) +
  geom_sf(data = st_cast(world, "MULTILINESTRING"),  size=0.1)+
  geom_sf_text(data = era5_4326_sf_pol, aes(label = DN), size=2) +
  coord_sf(xlim = c(-67.758632089, -66.741367911), ylim = c(-4.671851259, -3.828148741), expand = FALSE) +
  xlab("") + 
  ylab("") + 
  ggtitle("") +
  theme(panel.background = element_rect(fill = "aliceblue")) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
  theme(axis.ticks = element_blank()) +
  theme(plot.margin=grid::unit(c(0,0,0,0.2),"cm")) + 
  theme(panel.border = element_rect(colour = "red"))+
  theme(axis.ticks.length=unit(0, "null")) +
  theme(axis.ticks.margin=unit(0, "null")) +
  theme(axis.title.x=element_blank()) +
  theme(axis.title.y=element_blank()) +
  theme(plot.title = element_blank())

grid.arrange(big, arrangeGrob(corner1lt, corner2rt, corner3lb, corner4rb), ncol=2, widths=c(2.2,1))
```


```{r, include=FALSE}
#Get time serie from station 1 -Using ncdf4
# t.units <- ncatt_get(ncin, "time", "units")
# t.units
# time.array <-  ncvar_get(ncin, "time")
# nt = dim(time.array)
# nt
# tustr <- strsplit(t.units$value, " ")
# 
# #Our data is in Gregorian Calendar: "hours"      "since"      "1900-01-01" "00:00:00.0"
# anio_mes_dia = unlist(tustr)[3]
# 
# #Note that the dates are converted first to seconds!
# timestamp = as_datetime(c(time.array*60*60),origin=anio_mes_dia)
# timestamp_string <- as.character(timestamp)
# 
# statera5_xts = na.omit(xts(x= ncvar_get(ncin, 'fg10', start=c(1,1,1), count=c(1,1,ntime)), order.by = timestamp))
# 
# colnames(statera5_xts) = paste0("ERA5_", "1")

#Get time serie from station 1 -Using stars

time.array = read_ncdf(filename, var = c("time"))$time
t.units <- as.character(units(time.array)$numerator)
time.array = as.matrix(time.array)
nt = dim(time.array)
tustr <- strsplit(t.units, " ")
#Our data is in Gregorian Calendar: "hours"      "since"      "1900-01-01" "00:00:00.0"
anio_mes_dia = unlist(tustr)[3]
#Note that the dates are converted first to seconds!
timestamp = as_datetime(c(time.array*60*60),origin=anio_mes_dia)
timestamp_string <- as.character(timestamp)
fg10_1 = read_ncdf(filename, var = c(variablename), ncsub = cbind(start = c(1,1,1), count = c(1,1,ntime)))
statera5_xts = na.omit(xts(as.matrix(fg10_1$fg10), order.by = timestamp))

colnames(statera5_xts) = paste0("ERA5_", "1")
```

```{r, include=FALSE}
#calculate the taw periodogram from statera5_xts
plot(spectrum(statera5_xts))

#Calculate dayly mean and plot to so see if there is seasonal trend
statera5 = as.data.frame(statera5_xts)
#statera5$mytime = rownames(statera5)
statera5$mytime = timestamp
statera5$jday = as.numeric(format(statera5$mytime, '%j'))
Jday = 1:366
wind = as.matrix(statera5[1])
daymeans = sapply(split(wind, statera5$jday), mean)
plot(daymeans ~ Jday)
lines(lowess(daymeans ~ Jday, f = 0.1))

# subtract the trend:
meanwind = lowess(daymeans ~ Jday, f = 0.1)$y[statera5$jday]
velocity = apply(wind, 2, function(x) { x - meanwind })

#stl - Seasonal Decomposition of Time Series by Loess
#stl(statera5_xts, "per")
statera5.ts = ts(statera5[,1], frequency = 8760)
spectrum(statera5.ts, spans = c(3,5))
```

## Data Download and Data Organization

All data sources had different mechanisms for downloading. For IDEAM, the official procedure is through the e-mail _atencionalciudadano@ideam.gov.co_. For ISD all files are available in the FTP site ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/, organized in folders by years with _gzip_ files inside; there are many files by station (one file for each year available), with names in the format _ID-99999-YYYY.gz_, where ID is the USAF-ISD station identifier, and YYYY is the year. ERA5 data request uses a Python scripts with data size limit for download. Files with all IDEAM and ISD stations are available in Annex \@ref(rcode). For the Python code and commands to join NetCDF files of ERA5 data source see the Annex \@ref(datadownload). For data organization see the Annex \@ref(dbstoring).

\clearpage