Acceso_Datos_ERDDAP.Rmd

---
title: "Accediendo a datos de temperatura superficial desde ERDDAP"
author: "Denisse Fierro Arcos"
date: "2022-08-20"
output: 
  html_document: default
  github_document: default
---

Este script está basado en el tutorial desarrollado por Johnathan Evanilla para el OceanHackWeek 2022.


```{r}
suppressPackageStartupMessages({
  # from CRAN
  library(dplyr)
  library(stringr)
  library(raster)
  library(rerddap)
  library(ncdf4)
  library(ggplot2)
  library(mapdata)
  library(sf)
  library(leaflet)
  })

```



## Datos satelitales desde ERDDAP

[rerddap](https://docs.ropensci.org/rerddap/index.html)

Podemos obtener un listado de los datos disponibles a través de ERDAPP.

```{r}
s <- servers()
s
```
### Que tipo de datos esta disponibles en los servidores?
[Upwell ERDDAP Webpage](https://upwell.pfeg.noaa.gov/erddap/)

```{r}
ds_t <- ed_datasets(url="https://upwell.pfeg.noaa.gov/erddap/")
ds_g <- ed_datasets(url="https://upwell.pfeg.noaa.gov/erddap/", which="griddap")

paste("We have", sum(ds_t$tabledap!=""), "tabledap", sum(ds_g$griddap!=""), "griddap", sep=" ")
```
### Podemos hacer una busqueda para nuestra area de interes
Haremos un mapa con la extension de nuestra area de interes.

[leaflet](https://rstudio.github.io/leaflet/)

```{r}
lon <- c(-92, -89)
lat <- c(-2, 2)

leaflet() |>
  addTiles() |>
  setView(lng=sum(lon)/2, 
          lat=sum(lat)/2, 
          zoom=6) |>
  addRectangles(lng1 = lon[1],
                lng2 = lon[2],
                lat1 = lat[1],
                lat2 = lat[2])

```

## Busqueda avanzada de datos
Podemos poner limites espaciales y temporales

```{r}
out_td <- ed_search_adv(which = "griddap",
                        url="https://upwell.pfeg.noaa.gov/erddap/",
                        minLat=-2,
                        maxLat=2,
                        minLon=-92,
                        maxLon=-89,
                        minTime = "2004-01-01T00:00:00Z",
                        maxTime = "2022-07-31T00:00:00Z")

out_td
```

```{r}
paste("Found", length(out_td$alldata), "Datasets:", sep=" ")
```



### Buscando datos de temperatura de la superficie

```{r}
out <- ed_search(query = "SST", 
                 which = "griddap",
                 url="https://upwell.pfeg.noaa.gov/erddap/")

dataset_ids <- out$info$dataset_id

out
```
#Podemos hacer un mejor filtrado utilizando el tidyverse
Abajo ponemos un ejemplo de busqueda de datos MODIS obtenidos durante el dia a una resolucion temporal de un mes y resolucion espacial de 4 km.

```{r}
#Dataframe con informacion de datasets disponibles sobre temperatura superficial
sst_dataset <- out$info %>% 
  #Filtrado para datos MODIS
  filter(str_detect(str_to_lower(title), "modis") & 
           str_detect(str_to_lower(title), "daytime") &
           str_detect(str_to_lower(title), "4km") & 
           str_detect(str_to_lower(title), "monthly"))

sst_dataset
```


### Podemos obtener informacion adicional de los datos de nuestro interes

```{r}
info(sst_dataset$dataset_id[2], 
     url="https://upwell.pfeg.noaa.gov/erddap/")
```

### Accediendo datos desde el servidor

```{r}
#Podemos mantener los datos en la memoria solamente
temp_Gps <- griddap(x = sst_dataset$dataset_id[2],
                    latitude=c(-2, 2),
                    longitude=c(-92,-89),
                    time=c("2022-01-01","2022-02-28"),
                    fields = "sst",
                    #En este caso bajamos los datos como tipo csv
                    fmt = "csv")

#Podemos hacer cualquier cambio y guardarlos a nuestro equipo local
temp_Gps %>% 
  write.csv("Data/SST_Gps_2022-01-01_2022-02-28.csv", 
            row.names = F)
```


### Bajando datos desde el servidor

```{r}
#Podemos guardar datos directamente al disco utilizando la opcion 'store'
temp_Gps_nc <- griddap(x = sst_dataset$dataset_id[2],
                    latitude=c(-2, 2),
                    longitude=c(-92,-89),
                    time=c("2022-01-01","2022-02-28"),
                    fields = "sst",
                    #En este caso los datos seran de tipo netcdf
                    fmt = "nc",
                    #Podemos especificar la ruta donde guardremos los datos
                    store = disk(path = "Data/"))

temp_Gps_nc
```
Podemos ahora utilizar el paquete `raster` para subir los datos que acabamos de guardar como rasters.

```{r}
sst_raster <- raster("Data/e11136594a145f6c452c0e63764b231c.nc")
raster::plot(sst_raster)
```
Podemos usar el paquete `raster` para hacer calculos sencillos, como el promedio. 

```{r}
prom_sst_gps <- raster::mean(sst_raster, 3)

plot(prom_sst_gps)
```
Podemos tambien verificar el tipo de datos que hemos accedido.

```{r}
class(temp_Gps)
```
Asi como producir un resumen

```{r}
summary(temp_Gps)
```
Usando el tidyverse podemos hacer filtrado de datos. Por ejemplo, los datos de los meses de marzo y abril.

```{r}
temp_Gps %>% 
  filter(lubridate::month(time) == 3 |
           lubridate::month(time) == 4)
```

Graficando con ggplot

```{r}
temp_Gps %>% 
  filter(time == lubridate::ymd("2022-03-01")) %>% 
  as_tibble() %>% 
  ggplot(aes(x = longitude, y = latitude)) +
  geom_contour_filled(aes(z = sst)) +
  geom_sf(data = tierra, inherit.aes = F)+
  lims(x = c(-92, -89), y = c(-2, 2))
  
```