---
title: "Informe de copa del arbolado (HBC y dos paraboloides)"
author: "Felipe Bravo y Sergio Moro-Andrés"
date: "`r format(Sys.Date())`"
output: 
    html_document: 
    toc: true
    number_sections: true
---
# 1. Carga de paquetes de R

Solo usamos, además de R base, los paquetes dplyr y ggplot2 que ya están cargados en RStudio,
en el caso de que no lo estuvieran hay que correr las instrucciones intall.packages quitando la almohadilla

```{r}

#Si hay que instalar las librerias/paquetes se debe quitar la # de las siguientes 2 instrucciones
#install.packages(dplyr)
#install.packages(ggplot2)

library(dplyr)
library(ggplot2)

# Opciones base
options(dplyr.summarise.inform = FALSE)
```

# 1. Datos y preparación

Leemos el CSV en formato europeo (separador `;`, decimal `,`) con la instrucción `read.csv2` de  **base R**.
No usamos paquetes de limpieza de nombres; trabajamos con los nombres originales tal como están en el archivo.
El archivo Rmd y el archivo de datos deben estar en el mismo directorio. Por eso establecemos el directorio de trabajo ('C:/datosR'); 
cambia este directorio al que uses en tu ordenador y acuerdate cambiar el sentido de la barra, debe estar así '/' 
y no como viene por defecto en Windows (si usas Linux no tienes que preocuparte de esto)


```{r}
# establecemos el directorio de trabajo
setwd('C:/datosR')
getwd()

# Lectura usando base R: read.csv2 asume sep=';' y dec=','
df <- read.csv2("Treemetrics-total.csv", stringsAsFactors = FALSE, check.names = FALSE)

# Conversión de columnas numéricas relevantes (si estuvieran en texto)
num_cols <- c("tree_number","tree_height_field_m","dbh1_field_cm","dbh2_field_cm","dbh_field_cm",
              "tree_height_lidar_m","dbh_lidar_cm","CBH_lidar_m","Crown_length_lidar_m",
              "crown_ratio_lidar_percentage","Height_toMaximumCrownWidth_m",
              "Crown_widthX_lidar_m","Crown_widthY_lidar_m","Crown_width_lidar_m",
              "Crown_volume_lidarPlugin_m3","Crown_points_lidar_numbers","Crown_density_lidarPlugin_numbersbym3")

conv_num <- function(x){
  if(is.numeric(x)) return(x)
  as.numeric(gsub(",", ".", x))
}

for (cn in intersect(num_cols, names(df))) {
  df[[cn]] <- conv_num(df[[cn]])
}

# Normalizamos 'Crown Shape' a minúsculas
df[["Crown Shape"]] <- tolower(df[["Crown Shape"]])

# Vista previa
head(df, 5)
```

## 1.1 Reglas de omisión y coherencia

Comprobamos que los datos cumplen los requisitos necesarios para poder aplicar las funciones que vamos a definir despues: \([Volúmenes\ y\ áreas\ laterales\ de\ Copa]\) 
- Para el cálculo de copa se requieren (por árbol) que la siguientes variables cumplan los requisitos siguientes:
`Crown_widthX_lidar_m` (a, diámetro X) > 0, `Crown_widthY_lidar_m` (b, diámetro Y) > 0, `tree_height_lidar_m` (altura total) > 0, `CBH_lidar_m` (HBC) ≥ 0, `Height_toMaximumCrownWidth_m` (HmaxCW) ≥ 0.
- Para paraboloides, `HmaxCW` debe estar en el intervalo \([HBC,\ \text{altura total}]\).
- Si falta algo o no se cumple, **no se calcula** y se deja `NA`.


```{r}
has_crown_data <- function(row) {
  a      <- row[["Crown_widthX_lidar_m"]]
  b      <- row[["Crown_widthY_lidar_m"]]
  h_tot  <- row[["tree_height_lidar_m"]]
  hbc    <- row[["CBH_lidar_m"]]
  hmaxcw <- row[["Height_toMaximumCrownWidth_m"]]
  if (any(is.na(c(a, b, h_tot, hbc, hmaxcw)))) return(FALSE)
  (a > 0) && (b > 0) && (h_tot > 0) && (hbc >= 0) && (hmaxcw >= 0)
}
```

# 2. MÓDULO: **VOLÚMENES** (solo funciones)

Definiciones (con **diámetros** \(a,b\)):

- Área elíptica: \( A = \pi \frac{a}{2}\frac{b}{2} = \frac{\pi a b}{4} \).
- **Cono elíptico**: \( V = \frac{1}{3} A h = \frac{\pi a b h}{12} \), usando \( h_{\text{copa}} = \text{altura total} - \text{HBC} \).
- **Paraboloide elíptico** (un tramo): \( V = \frac{1}{2} A h = \frac{\pi a b h}{8} \).
- **Paraboloide en dos tramos**: \( V_{\text{total}} = V(h_{\text{bot}}) + V(h_{\text{top}}) \) con  
  \( h_{\text{bot}} = \text{HmaxCW} - \text{HBC} \),  
  \( h_{\text{top}} = \text{Altura total} - \text{HmaxCW} \).

```{r}
area_elipse <- function(a, b) {
  if (any(is.na(c(a, b))) || a <= 0 || b <= 0) return(NA_real_)
  pi * a * b / 4
}

vol_cone_eliptico <- function(a, b, altura_total, hbc) {
  A <- area_elipse(a, b)
  h_copa <- altura_total - hbc
  if (is.na(A) || is.na(h_copa) || h_copa <= 0) return(NA_real_)
  (1/3) * A * h_copa  # = pi*a*b*h_copa/12
}

vol_paraboloide_tramo <- function(a, b, h_tramo) {
  A <- area_elipse(a, b)
  if (is.na(A) || is.na(h_tramo) || h_tramo <= 0) return(NA_real_)
  0.5 * A * h_tramo   # = pi*a*b*h_tramo/8
}

vol_paraboloide_dos_tramos <- function(a, b, altura_total, hbc, hmaxcw) {
  if (any(is.na(c(a, b, altura_total, hbc, hmaxcw)))) return(NA_real_)
  h_bot <- hmaxcw - hbc
  h_top <- altura_total - hmaxcw
  if (a <= 0 || b <= 0 || h_bot <= 0 || h_top <= 0) return(NA_real_)
  vol_paraboloide_tramo(a, b, h_bot) + vol_paraboloide_tramo(a, b, h_top)
}
```

# 3. MÓDULO: **ÁREA LATERAL** (solo funciones)

- **Radio equivalente por área** (círculo con igual área que la elipse): \( r_{eq} = \sqrt{ab}/2 \) (con **diámetros** \(a,b\)).  
- **Semiparaboloide** (apex en \(z=0\), altura \(h\)): \( r(z) = r_{eq}\sqrt{z/h} \).  
- Área lateral por **frustas**: discretización de \([0,h]\) en \(n\) segmentos; para cada segmento  
  generatriz \( s = \sqrt{(\Delta z)^2 + (\Delta r)^2} \), radio medio \( r_m \),  
  \( A_{\text{seg}} \approx 2\pi r_m s \); suma total.  
- **Cono**: \( SA \approx \pi r_{eq} l \), con \( l = \sqrt{r_{eq}^2 + h_{\text{copa}}^2} \).

```{r}
r_eq_from_ab <- function(a, b) {
  if (any(is.na(c(a, b))) || a <= 0 || b <= 0) return(NA_real_)
  sqrt(a * b) / 2
}

sa_paraboloid_part <- function(r_eq, h, n = 200) {
  if (is.na(r_eq) || is.na(h) || r_eq <= 0 || h <= 0) return(NA_real_)
  z <- seq(0, h, length.out = n + 1)
  r <- r_eq * sqrt(z / h)
  dz <- diff(z)
  dr <- diff(r)
  s  <- sqrt(dz^2 + dr^2)
  r_m <- (r[-1] + r[-length(r)]) / 2
  sum(2 * pi * r_m * s)
}

sa_paraboloid_two_parts <- function(a, b, altura_total, hbc, hmaxcw, n = 200) {
  if (any(is.na(c(a, b, altura_total, hbc, hmaxcw)))) return(NA_real_)
  r_eq <- r_eq_from_ab(a, b)
  h_bot <- hmaxcw - hbc
  h_top <- altura_total - hmaxcw
  if (is.na(r_eq) || r_eq <= 0 || h_bot <= 0 || h_top <= 0) return(NA_real_)
  sa_paraboloid_part(r_eq, h_bot, n) + sa_paraboloid_part(r_eq, h_top, n)
}

sa_cone_eliptico <- function(a, b, altura_total, hbc) {
  r_eq <- r_eq_from_ab(a, b)
  h_copa <- altura_total - hbc
  if (is.na(r_eq) || is.na(h_copa) || r_eq <= 0 || h_copa <= 0) return(NA_real_)
  l <- sqrt(r_eq^2 + h_copa^2)
  pi * r_eq * l
}
```

# 4. Cálculo por árbol e integración

Aplicamos las funciones por árbol con **dplyr**. Mantendremos `NA` si faltan datos o si `HmaxCW` no está dentro del rango \([HBC,\ \text{altura total}]\).

```{r}
calc <- df %>%
  rowwise() %>%
  mutate(
    # Accesos cortos
    a = Crown_widthX_lidar_m,
    b = Crown_widthY_lidar_m,
    h_tot  = tree_height_lidar_m,
    hbc    = CBH_lidar_m,
    hmaxcw = Height_toMaximumCrownWidth_m,
    shape  = tolower(`Crown Shape`),

    # Validaciones
    valid_data = has_crown_data(cur_data()),
    valid_hmax = ifelse(is.na(hmaxcw) | is.na(hbc) | is.na(h_tot), FALSE,
                        (hmaxcw >= hbc) & (hmaxcw <= h_tot)),

    # Alturas según tus definiciones
    h_copa = ifelse(valid_data, h_tot - hbc, NA_real_),
    h_bot  = ifelse(valid_data & valid_hmax, hmaxcw - hbc, NA_real_),
    h_top  = ifelse(valid_data & valid_hmax, h_tot - hmaxcw, NA_real_),

    # Volumen
    volume_geom_m3 = case_when(
      valid_data & shape == "paraboloid" & valid_hmax ~ vol_paraboloide_dos_tramos(a, b, h_tot, hbc, hmaxcw),
      valid_data & shape == "conical"                 ~ vol_cone_eliptico(a, b, h_tot, hbc),
      TRUE                                            ~ NA_real_
    ),

    # Área lateral
    lateral_area_geom_m2 = case_when(
      valid_data & shape == "paraboloid" & valid_hmax ~ sa_paraboloid_two_parts(a, b, h_tot, hbc, hmaxcw, n = 200),
      valid_data & shape == "conical"                 ~ sa_cone_eliptico(a, b, h_tot, hbc),
      TRUE                                            ~ NA_real_
    ),

    # Radio equivalente (referencia)
    r_eq_max_m = r_eq_from_ab(a, b)
  ) %>%
  ungroup()

# Vista previa de resultados
head(calc[, c("Species","Crown Shape","h_copa","h_bot","h_top","volume_geom_m3","lateral_area_geom_m2")], 10)
```

# 5. Resúmenes

```{r}
# Resumen global vs volumen del plugin (si existe)
resumen_global <- calc %>%
  summarise(
    n = n(),
    vol_plugin_mean = mean(Crown_volume_lidarPlugin_m3, na.rm = TRUE),
    vol_geom_mean   = mean(volume_geom_m3, na.rm = TRUE),
    bias_vol        = mean(volume_geom_m3 - Crown_volume_lidarPlugin_m3, na.rm = TRUE),
    rmse_vol        = sqrt(mean((volume_geom_m3 - Crown_volume_lidarPlugin_m3)^2, na.rm = TRUE))
  )

resumen_especie <- calc %>%
  group_by(Species) %>%
  summarise(
    n = n(),
    h_copa_mean = mean(h_copa, na.rm = TRUE),
    h_bot_mean  = mean(h_bot, na.rm = TRUE),
    h_top_mean  = mean(h_top, na.rm = TRUE),
    a_mean      = mean(Crown_widthX_lidar_m, na.rm = TRUE),
    b_mean      = mean(Crown_widthY_lidar_m, na.rm = TRUE),
    r_eq_mean   = mean(r_eq_max_m, na.rm = TRUE),
    vol_geom_mean = mean(volume_geom_m3, na.rm = TRUE),
    sa_geom_mean  = mean(lateral_area_geom_m2, na.rm = TRUE)
  )

# Impresión base
resumen_global
resumen_especie
```

# 6. Figuras y exportación de PNG

```{r}
# Volumen: plugin vs geométrico
p_vol <- ggplot(calc, aes(x = Crown_volume_lidarPlugin_m3, y = volume_geom_m3, color = Species)) +
  geom_point(size = 3) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  labs(x = "Volumen plugin LIDAR (m³)", y = "Volumen geométrico (m³)",
       title = "Volumen de copa: plugin vs geométrico") +
  theme_minimal()

# Área lateral vs longitud total de copa (h_copa)
p_sa <- ggplot(calc, aes(x = h_copa, y = lateral_area_geom_m2, color = `Crown Shape`)) +
  geom_point(size = 3) +
  labs(x = "Longitud total de copa h_copa (m)", y = "Área lateral (m²)",
       title = "Área lateral de copa vs longitud de copa") +
  theme_minimal()

# Altura vs diámetro

h_d <- ggplot(calc, aes(x = dbh_field_cm, y = tree_height_field_m, color = `Crown Shape`)) +
  geom_point(size = 3) +
  labs(x = "diámetro normal (cm)", y = "Áltura total (m)",
       title = "Relación altura-diámetro (mediciones tradicionales") +
  theme_minimal()

# Diámetro forcípula vs diámetro LiDAR

dtrad_dlidar <- ggplot(calc, aes(x = dbh_field_cm, y = dbh_lidar_cm, color = `Crown Shape`)) +
  geom_point(size = 3) +
  labs(x = "Diámetro normal (cm) con forcípula", y = "Diámetro normal (cm) con LiDAR",
       title = "Relación entre diámetros con forcípula y LiDAR") +
  theme_minimal()

# Altura tradicional vs Altura LiDAR

htrad_hlidar <- ggplot(calc, aes(x = tree_height_field_m, y = tree_height_lidar_m, color = `Crown Shape`)) +
  geom_point(size = 3) +
  labs(x = "Altura total (m) con hipsómetro", y = "Áltura total (m) con LiDAR",
       title = "Relación entre alturas con hipsómetro y LiDAR") +
  theme_minimal()

p_vol
p_sa
h_d
dtrad_dlidar
htrad_hlidar

# Guardar PNG
if (!dir.exists("figuras")) dir.create("figuras") # si no existe el directorio 'figuras' lo creamos
ggsave("figuras/volumen_plugin_vs_geom.png", p_vol, width = 8, height = 6, dpi = 300)
ggsave("figuras/area_lateral_vs_longitud_copa.png", p_sa, width = 8, height = 6, dpi = 300)
ggsave("figuras/relacion_altura_diámetro.png", h_d, width = 8, height = 6, dpi = 300)
ggsave("figuras/relacion_diametro-trad_diametro-lidar.png", dtrad_dlidar, width = 8, height = 6, dpi = 300)
ggsave("figuras/relacion_altura-trad_altura-lidar.png", htrad_hlidar, width = 8, height = 6, dpi = 300)

# Mostrar rutas de imágenes guardadas (impresión base)
list.files("figuras", pattern = "\\.png$", full.names = TRUE)
```

# 7. Diagnóstico

```{r}
# Árboles sin cálculo (NA) o HmaxCW fuera de rango
diag <- calc %>%
  mutate(id = seq_len(n())) %>%
  filter(is.na(volume_geom_m3) | is.na(lateral_area_geom_m2))

# Impresión base
subset(diag, select = c("id","Species","Crown Shape","tree_height_lidar_m","CBH_lidar_m","Height_toMaximumCrownWidth_m","h_copa","h_bot","h_top"))
```

# 8. Reproducibilidad

```{r}
sessionInfo()
```
