convertir_raw2nc.qmd

---
title: "Convertir archivos raw (EK80) a netCDF"
format: html
editor: visual
---

## 1. Primero convertir raw -\> nc con python

Instalar paquetes en terminal con `pip install 'package'`. NOTA: `cartopy` no se pudo instalar. Una vez instalados, cargar en la sesión.

```{python}
from pathlib import Path

import fsspec
import numpy as np
import geopandas as gpd
import xarray as xr

import matplotlib.pyplot as plt
from shapely.geometry import box
#import cartopy.crs as ccrs
#import cartopy.io.img_tiles as cimgt
#from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import echopype as ep
from echopype.qc import exist_reversed_time

import warnings
warnings.simplefilter("ignore", category=DeprecationWarning)
```

Definir directorio con archivos raw

```{python}
rawdirpath = "raw"
```

Nombres de archivos raw

```{python}
import glob
s3rawfiles = glob.glob(f"{rawdirpath}/*.raw")
s3rawfiles
```

```{python}
# Parámetros ambientales
env_params = {
    'temperature': 26.15,   # temperature in degree Celsius
    'salinity': 34.96,     # salinity in PSU
    'pressure': 25     # pressure in dbar
}
```

Abrir raw como como EchoData y convertir a netCDF

```{python}
for i in s3rawfiles:
  ed = ep.open_raw(i, sonar_model='EK80')
  ed.to_netcdf(save_path='./converted')
  ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode = "CW", encode_mode = "complex")
  # Reduce data based on sample number
  ds_MVBS = ep.preprocess.compute_MVBS_index_binning(
     ds_Sv,             # calibrated Sv dataset
     range_sample_num=30,  # number of sample bins to average along the range_sample dimensionm
     ping_num=5         # number of pings to average
  )
  ii = i.replace("raw\\", "")
  ii = ii.replace(".raw", "")
  nombre = "./converted/"+str(ii) + "_MVBS" + ".nc"
  #print(nombre)
  ds_MVBS.to_netcdf(nombre)
```

## 2. Importar netCDF simplificado en R

```{r}
library(ncdf4)
nc <- list.files(path = "./converted", pattern = glob2rx("*MVBS.nc"), full.names = TRUE)
```

Abrir un archivo

```{r}
ncf <- nc_open(nc[1])
```

Inspección del contenido

```{r}
names(ncf$var)
```

Extraer lo necesario para crear el objeto de clase echogram

```{r}
Sv <- ncvar_get(ncf, "Sv") # extraer valores de Sv
f1 <- Sv[, , 1] # solo 38 kHz (primera tabla de la 3a dimensión)
#f1 <- f1[nrow(f1):1, ] # ordenar matriz por orden descendente de filas
image(t(f1)) # primera visualización
```

Ahora `pingTime`

```{r}
pt <- ncf$dim$ping_time$vals
pt <- as.POSIXct(pt, tz = "UTC", format = "%Y-%m-%d %H:%M:%OS", origin = "1900-01-01  00:00:00")
```

Vector de profundidades de las muestras

```{r}
# calcular manualmente sample range (profundidad de las muestras)
# número de muestras
nd <- nrow(f1)
# velocidad del sonido
ss <- 1537.93 # buscarla en los datos calibrados
# sample interval
si <- 0.000512
# calcular sample length
sl <- (ss * si) / 2
R <- rep(sl, nd) # it must be the same for all pings
R <- as.numeric(cumsum(R)) # sample range  
```

Crear objeto clase echogram

```{r}
attr(f1, "frequency") <- "38 kHz"
eco <- list(depth = R,
            Sv = f1,
            pings = data.frame(
              pingTime = pt, # clase POSIXct
              detBottom = NA,
              speed = NA,
              cumdist = NA
            ))
class(eco) <- "echogram"
```

Figura

```{r}
echogram(eco)
```

## 3. Código R (archivos netCDF crudos)

Cargar paquete ncdf4

```{r}
library(ncdf4)
nc <- list.files(path = "./converted", pattern = glob2rx("*.nc"), full.names = TRUE)
```

Abrir un archivo

```{r}
ncf <- nc_open(nc[1])
```

Inspección del contenido

```{r}
names(ncf$var)

# verificar que existen
#ncvar_get(ncf, "Environment/absorption_indicative")
ncvar_get(ncf, "Platform/frequency_nominal")
lat <- ncvar_get(ncf, "Platform/latitude")
lon <- ncvar_get(ncf, "Platform/longitude")

plot(lon, lat, asp = 1)
```

Matriz de datos

```{r}
Sv <- ncvar_get(ncf, "Sonar/Beam_group1/backscatter_r")
dim(Sv)
f1 <- Sv[1, , , 1]
f1 <- f1[nrow(f1):1, ]
image(t(f1))
```

## Uso de funcion `getPingTimeEK80()`

```{r}
source("getPingTimeEK80.R")
pt <- getPingTimeEK80(nc[2])
pt[1:10]
# ver fracciones de segundo 
strftime(pt[1:10], "%Y-%m-%d %H:%M:%OS2")
```

Ecuación para convertir *Received power* (*Pr*) en $S_v$ (EK60)

$$
S_V(R, P_r) = P_r + 20log(R) + 2\alpha R - 10log(\frac{P_t G_0^2 \lambda ^2}{16 \pi^2}) -10log(\frac{c \tau \psi}{2}) -2 S_a\;corr
$$

Ec. para convertir *Pr* en $TS$

$$
TS(R, P_r) = P_r + 40log(R) + 2\alpha R - 10log(\frac{P_t G_0^2 \lambda ^2}{16 \pi^2})
$$

## Convertir a echogram

instalar el paquete

```{r}
# OJO SOLO UNA VEZ. Despues comentar las dos lineas con #
library(devtools)
install_github("hvillalo/echogram")
```

Ver la estructura de la clase

```{r}
hacfile <- system.file("extdata", "D20150510-T202221.hac", package = "echogram")
echo1 <- read.echogram(hacfile, channel = 1)
str(echo1)
```

Crear obj echogram desde el nc

```{r}
eco <- list(depth = ,
            Sv = , # matriz
            pings = data.frame(
              pingTime = , # clase POSIXct
              detBottom = NA,
              speed = NA,
              cumdist = NA
            ))
```