BE Group – Recruitment Test Report

Este relatório foi desenvolvido como parte do processo seletivo da BE Group, conforme instruções disponibilizadas no repositório oficial:

🔗 https://github.com/pecard/testeR

Contexto

O teste consiste em dois desafios analíticos envolvendo:

  • Análise espacial de produção de energia eólica em Portugal continental;
  • Análise de movimentos de aves em relação à maré e período do dia.

O objetivo principal é demonstrar:

  • Capacidade de estruturação analítica;
  • Organização de workflow em R;
  • Tratamento e integração de dados espaciais e temporais;
  • Aplicação de estatística descritiva e inferencial;
  • Produção de visualizações claras e interpretáveis.

1. Teste 1 - Wind Energy Production in Portugal

Introdução

Este relatório apresenta a análise espacial e temporal das instalações de energia eólica em Portugal continental, com foco em responder o Teste 1 de avaliação de dados e visualização geoespacial.
O objetivo principal é explorar, analisar e visualizar a distribuição das usinas eólicas considerando diversos fatores ambientais e administrativos.

Perguntas do Teste 1

As análises abaixo foram desenvolvidas para responder às seguintes questões:

  1. Mapear as instalações eólicas em Portugal continental considerando os limites do país e dos distritos. Aplicar legenda apropriada e escolher um nível de agregação para as usinas.
  2. Mapear as instalações eólicas e sua distribuição em relação à altitude e áreas protegidas.
  3. Selecionar o tema adequado para o Modelo Digital de Elevação (DEM).
  4. Determinar a capacidade instalada nacional total (em GW) em 2019.
  5. Avaliar como as instalações estão distribuídas em relação à altitude, plotando a distribuição altitudinal e descrevendo os resultados.
  6. Analisar o crescimento da capacidade instalada (em GW) desde a primeira instalação.
  7. Verificar se as áreas protegidas são afetadas, identificando quais e em que proporção (número de sítios Natura 2000).

1.1. Carregamento dos Dados e Pacotes

# ===============================
# Carregamento dos pacotes e funções
# ===============================
library(tidyverse)  # Manipulação de dados
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.3
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   4.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)         # Spatial features
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(terra)      # Raster handling
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.5
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(ggspatial)  # Map layers e escala
## Warning: package 'ggspatial' was built under R version 4.3.3
library(leaflet)    # Mapas interativos
## Warning: package 'leaflet' was built under R version 4.3.3
library(viridis)    # Paletas de cores
## Warning: package 'viridis' was built under R version 4.3.3
## Carregando pacotes exigidos: viridisLite
## Warning: package 'viridisLite' was built under R version 4.3.3
library(ggnewscale) # Múltiplas escalas de cores no ggplot

# Funções customizadas
source("R/01_setup.R")
## Warning: package 'data.table' was built under R version 4.3.3
## 
## Attaching package: 'data.table'
## 
## The following object is masked from 'package:terra':
## 
##     shift
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## Warning: package 'raster' was built under R version 4.3.3
## Carregando pacotes exigidos: sp
## Warning: package 'sp' was built under R version 4.3.3
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:viridis':
## 
##     viridis_pal
## 
## The following object is masked from 'package:terra':
## 
##     rescale
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## Warning: package 'car' was built under R version 4.3.3
## Carregando pacotes exigidos: carData
## Warning: package 'carData' was built under R version 4.3.3
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## Carregando pacotes exigidos: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:raster':
## 
##     getData
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## This is mgcv 1.9-0. For overview type 'help("mgcv-package")'.
## 
## Attaching package: 'MASS'
## 
## The following objects are masked from 'package:raster':
## 
##     area, select
## 
## The following object is masked from 'package:terra':
## 
##     area
## 
## The following object is masked from 'package:dplyr':
## 
##     select
source("R/functions.R")

# ===============================
# Carregamento dos dados
# ===============================
load("data/e2p.rda")      # Instalações eólicas
load("data/adm0_pt.rda")  # Limite país
load("data/adm1_pt.rda")  # Limite distritos
load("data/rn2000.rda")   # Áreas Natura 2000

# ===============================
# Opções globais dos chunks
# ===============================
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)

1.2. Limpeza, ETL e Padronização dos Dados

# ===============================
# Limpeza e preparação dos dados
# ===============================

# Renomear colunas e converter tipos
e2p_clean <- e2p %>%
  rename(
    latitude = latitude_wgs84_o,
    longitude = longitude_wgs84_o,
    potencia_mw = potencia_instalada_mw,
    ano = ano_de_entrada_em_funcionamento
  ) %>%
  mutate(
    across(c(latitude, longitude, potencia_mw), ~ as.numeric(gsub(",", ".", .))),
    ano = as.numeric(ano)
  )

# ===============================
# Converter para sf e definir CRS PT-TM06
# ===============================
e2p_sf <- st_as_sf(e2p_clean, coords=c("longitude", "latitude"), crs=4326)
adm0_pt <- st_transform(adm0_pt, 3763)
adm1_pt <- st_transform(adm1_pt, 3763)
e2p_sf <- st_transform(e2p_sf, 3763)

# ===============================
# Filtrar apenas instalações em terra
# ===============================
e2p_mainland <- e2p_sf[st_intersects(e2p_sf, adm0_pt, sparse=FALSE), ]

1.3. Distribuição das Instalações Eólicas em Portugal Continental

O mapa apresenta as instalações de energia eólica localizadas em Portugal continental, sobrepondo os limites do país e dos distritos. As instalações estão classificadas em quatro classes de capacidade instalada: <5 MW, 5–20 MW, 20–50 MW e >50 MW. Este agrupamento permite observar visualmente onde estão concentradas as maiores e menores usinas eólicas, facilitando análises de planejamento energético e impacto territorial.

# Criando classe por potência instalada dos empreendimentos (MW)
e2p_sf <- e2p_sf %>%
  mutate(
    classe_potencia = cut(
      potencia_mw,
      breaks = c(-Inf, 5, 20, 50, Inf),
      labels = c("< 5 MW", "5–20 MW", "20–50 MW", "> 50 MW"),
      right = FALSE
    )
  )

# Filtrar apenas instalações em terra
e2p_mainland <- e2p_sf[st_intersects(e2p_sf, adm0_pt, sparse=FALSE), ]

# Plotando o mapa
ggplot() +
  geom_sf(data = adm0_pt, fill = "grey95", color = "black") +
  geom_sf(data = adm1_pt, fill = NA, color = "grey70", size = 0.3) +
  geom_sf(data = e2p_mainland,
          aes(color = classe_potencia),
          size = 2, alpha = 0.8) +
  scale_color_viridis_d(name = "Capacidade Instalada (MW)") +
  labs(
    title = "Instalações Eólicas em Portugal Continental",
    subtitle = "Classificação das instalações por capacidade instalada",
    caption = "Projeção: PT-TM06 (EPSG:3763)\nFonte: Base de dados de instalações eólicas"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    style = north_arrow_fancy_orienteering
  ) +
  theme_minimal() +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold")
    )

1.4. Distribuição das Instalações Eólicas em Relação à Altitude e Áreas Protegidas

O mapa apresenta a distribuição das instalações eólicas em Portugal continental, considerando a altitude do terreno (DEM de 25 m reamostrado) e as áreas protegidas Natura 2000.

O relevo é representado pelo DEM (Digital Elevation Model) de 25m, com cores graduadas do verde claro ao azul escuro indicando a elevação.

As áreas protegidas são destacadas em verde claro sem contorno, permitindo observar quais regiões estão sob proteção ambiental.

As usinas eólicas são representadas por pontos pretos, indicando sua localização exata.

Esta visualização permite analisar a relação entre a instalação de parques eólicos e o relevo, além de identificar possíveis impactos em áreas de conservação. impactos em áreas de conservação.

# Extrair altitude do DEM 25m reamostrado
dem <- rast("data/dem_srtm_pt_25m.tif")
crs(dem) <- "EPSG:3763"

e2p_vect <- vect(e2p_mainland)
e2p_mainland$elevation_m <- extract(dem, e2p_vect)[,2]

# Preparar DEM para visualização mais leve
dem_vis <- aggregate(dem, fact = 8)
## |---------|---------|---------|---------|=========================================                                          
dem_vis <- mask(dem_vis, vect(adm0_pt))
dem_df <- as.data.frame(dem_vis, xy = TRUE, na.rm = TRUE)
colnames(dem_df)[3] <- "elevation"

# Garantir CRS correto para áreas protegidas
rn2000 <- st_transform(rn2000, 3763)

# Plot final
ggplot() +
  # DEM
  geom_raster(data = dem_df,
              aes(x = x, y = y, fill = elevation),
              alpha = 0.9) +
  scale_fill_viridis_c(name = "Altitude (m)", option = "C") +

  # Resetar escala fill
  new_scale_fill() +

  # Áreas Natura 2000
  geom_sf(data = rn2000,
          aes(fill = "Áreas Protegidas"),
          color = "darkgreen",
          alpha = 0.5) +
  scale_fill_manual(name = "", values = c("Áreas Protegidas" = "lightgreen")) +

  # Resetar escala color
  new_scale_color() +

  # Instalações eólicas
  geom_sf(data = e2p_mainland,
          aes(color = "Usinas Eólicas"),
          size = 2,
          alpha = 0.9) +
  scale_color_manual(name = "", values = c("Usinas Eólicas" = "black")) +

  # Limites administrativos
  geom_sf(data = adm0_pt, fill = NA, color = "black", linewidth = 0.6) +
  geom_sf(data = adm1_pt, fill = NA, color = "grey60", linewidth = 0.3) +

  labs(
    title = "Distribuição das Instalações Eólicas, Altitude e Áreas Protegidas",
    subtitle = "Visualização das usinas em relação ao relevo e aos sítios Natura 2000",
    caption = "DEM de 25 m reamostrado\nProjeção: PT-TM06 (EPSG:3763)\nFonte: Conjunto de dados de teste do BE Group"
  ) +
  annotation_scale(location = "bl", width_hint = 0.3) +
  annotation_north_arrow(
    location = "br",
    which_north = "true",
    style = north_arrow_fancy_orienteering
  ) +
  coord_sf(crs = 3763) +
  guides(
    color = guide_legend(order = 1),   # Usinas
    fill = guide_legend(order = 2),    # Áreas protegidas
    fill = guide_colorbar(order = 3)   # Altitude
  ) +
  theme_minimal()

1.5. Qual foi a capacidade instalada nacional total (em GW) em 2019?

Capacidade instalada nacional em 2019: A capacidade total de geração eólica em Portugal continental até o final de 2019 foi de aproximadamente 5.37 GW.

Distribuição altitudinal: As instalações eólicas estão localizadas em altitudes entre 2 e 1353 metros, com maior concentração em regiões próximas de 300–400 m e 900–1000 m, conforme indicado pelo histograma, mediana e quartis da distribuição. Isso reflete a preferência por áreas com ventos mais consistentes e relevo adequado à instalação de turbinas.

# ============================================
# Capacidade instalada nacional em 2019 (GW)
# ============================================

# Remover geometria para somar MW
dados_sem_geom <- st_drop_geometry(e2p_mainland)

total_2019 <- dados_sem_geom %>%
  filter(!is.na(ano), ano <= 2019) %>%
  summarise(total_mw = sum(potencia_mw, na.rm = TRUE)) %>%
  mutate(total_gw = total_mw / 1000)

cat(
  "Capacidade instalada nacional em 2019:",
  round(total_2019$total_gw, 2),
  "GW\n"
)
## Capacidade instalada nacional em 2019: 5.37 GW
# ============================================
# Distribuição altitudinal das instalações
# ============================================

ggplot(e2p_mainland, aes(x = elevation_m)) +
  geom_histogram(aes(y = ..density..),
                 bins = 30,
                 fill = "grey70",
                 color = "white") +
  geom_density(color = "red", linewidth = 1) +
  labs(
    title = "Distribuição altitudinal das instalações eólicas",
    x = "Altitude (m)",
    y = "Densidade",
    caption = "O histograma mostra a frequência das altitudes das usinas. A curva vermelha indica a densidade estimada."
  ) +
  theme_minimal()

# Estatísticas resumidas
summary(e2p_mainland$elevation_m)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     2.0   360.5   902.0   753.4  1065.5  1353.0

1.6. Como a capacidade instalada (em GW) evoluiu desde a primeira instalação?

O gráfico apresenta a evolução acumulada da capacidade instalada de energia eólica em Portugal continental desde a primeira usina, instalada em 1992, até 2019. Inicialmente, observa-se um crescimento lento, seguido por uma fase de forte expansão a partir do início dos anos 2000, refletindo o aumento do investimento em energia renovável. Por volta de 2015, a expansão tende a se estabilizar, mostrando um comportamento geral de curva sigmoidal. A linha tracejada indica o ano da primeira instalação, servindo como referência temporal para toda a evolução.

# Filtrar crescimento acumulado e remover NA's
growth <- e2p_mainland %>%
  filter(!is.na(ano)) %>%
  group_by(ano) %>%
  summarise(total_mw = sum(potencia_mw, na.rm = TRUE)) %>%
  arrange(ano) %>%
  mutate(
    cum_mw = cumsum(total_mw),
    cum_gw = cum_mw / 1000
  )

# Plot do crescimento acumulado
ggplot(growth, aes(x = ano, y = cum_gw)) +
  geom_line(linewidth = 1.1, color = "darkblue") +
  geom_point(size = 2, color = "darkblue") +
  
  # Indicar a primeira instalação
  geom_vline(xintercept = min(growth$ano),
             linetype = "dashed",
             color = "grey50") +
  annotate(
    "text",
    x = min(growth$ano) + 1,
    y = max(growth$cum_gw) * 0.1,
    label = "Primeira instalação",
    hjust = 0,
    size = 3
  ) +
  
  labs(
    title = "Crescimento Acumulado da Capacidade Instalada de Energia Eólica em Portugal Continental",
    subtitle = "Evolução desde a primeira instalação até 2019",
    x = "Ano",
    y = "Capacidade Instalada (GW)"
  ) +
  
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  )

1.7. As áreas protegidas são afetadas? Onde e em que extensão (número de sítios)?

O mapeamento das instalações eólicas em relação aos sítios Natura 2000 permite identificar as áreas de conservação potencialmente impactadas. No total, 70 usinas eólicas estão localizadas dentro de áreas protegidas, representando aproximadamente 29,8 % do total de usinas.

Em termos de sítios Natura 2000, 13 áreas apresentam pelo menos uma instalação eólica em seu interior, correspondendo a cerca de 21,3 % do total de sítios.

No gráfico a seguir:

  • Cinza claro: todas as áreas Natura 2000;
  • Vermelho: sítios Natura 2000 que contêm pelo menos uma usina;
  • Preto: usinas localizadas dentro de áreas protegidas.

Esta visualização facilita a avaliação espacial do impacto potencial sobre áreas de conservação e fornece informações importantes para orientar políticas de mitigação e estratégias de monitoramento ambiental.

# ===============================
# 1. Identificar usinas dentro de áreas protegidas
# ===============================

# Matriz de interseção (usinas vs Natura 2000)
intersections_wind <- st_intersects(e2p_mainland, rn2000)

# Criar indicador lógico
e2p_mainland$inside_protected <- lengths(intersections_wind) > 0

# Número de usinas dentro
n_wind_inside <- sum(e2p_mainland$inside_protected)

# Proporção
prop_wind_inside <- n_wind_inside / nrow(e2p_mainland) * 100

# ===============================
# 2. Identificar sítios Natura 2000 afetados
# ===============================

intersections_sites <- st_intersects(rn2000, e2p_mainland)

# Sítios que contêm pelo menos uma usina
sites_affected <- rn2000[lengths(intersections_sites) > 0, ]

# Número e proporção de sítios afetados
n_sites_affected <- nrow(sites_affected)
prop_sites_affected <- n_sites_affected / nrow(rn2000) * 100

# ===============================
# 3. Printar resultados
# ===============================
cat("Usinas dentro de áreas protegidas:", n_wind_inside, "\n")
## Usinas dentro de áreas protegidas: 70
cat("Percentual de usinas dentro de áreas protegidas:", round(prop_wind_inside, 2), "%\n\n")
## Percentual de usinas dentro de áreas protegidas: 29.79 %
cat("Sítios Natura 2000 afetados:", n_sites_affected, "\n")
## Sítios Natura 2000 afetados: 13
cat("Percentual de sítios afetados:", round(prop_sites_affected, 2), "%\n")
## Percentual de sítios afetados: 21.31 %
# ===============================
# 4. Plotar mapa
# ===============================

ggplot() +

  # Limite nacional
  geom_sf(data = adm0_pt,
          fill = "grey95",
          color = "black",
          linewidth = 0.6) +

  # Distritos
  geom_sf(data = adm1_pt,
          fill = NA,
          color = "grey70",
          linewidth = 0.3) +

  # Natura 2000 (todas)
  geom_sf(data = rn2000,
          fill = "grey85",
          color = NA,
          alpha = 0.5) +

  # Natura 2000 afetadas
  geom_sf(data = sites_affected,
          fill = "red",
          color = "darkred",
          alpha = 0.6) +

  # Usinas dentro de áreas protegidas
  geom_sf(data = e2p_mainland[e2p_mainland$inside_protected, ],
          color = "black",
          size = 2) +

  coord_sf(crs = 3763) +

  labs(
    title = "Sítios Natura 2000 Afetados por Usinas Eólicas",
    subtitle = "Usinas localizadas dentro de áreas protegidas",
    caption = "Projeção: PT-TM06 (EPSG:3763)"
  ) +

  theme_minimal()

Desafio 1 - Mapa Interativo das Instalações Eólicas

O mapa a seguir apresenta todas as instalações eólicas em Portugal continental de forma interativa. Cada ponto representa uma usina eólica, e ao clicar nele é possível visualizar informações detalhadas, incluindo o nome da instalação, o ano de entrada em operação e a capacidade instalada (em MW).

O tamanho dos pontos é proporcional à potência instalada, permitindo identificar rapidamente as usinas de maior capacidade. Essa visualização interativa facilita a exploração espacial dos dados e a análise de distribuição das instalações.

# Garantir WGS84
e2p_map <- st_transform(e2p_mainland, 4326)

# Criar popups
e2p_map <- e2p_map %>%
  mutate(
    popup_info = paste0(
      "<b>Instalação:</b> ", nome, "<br>",
      "<b>Ano:</b> ", ano, "<br>",
      "<b>Capacidade instalada:</b> ", round(potencia_mw, 1), " MW"
    )
  )

# Mapa interativo
leaflet(e2p_map) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    radius = ~sqrt(potencia_mw) / 2,   # escala visual pela potência
    color = "darkred",
    stroke = FALSE,
    fillOpacity = 0.8,
    popup = ~popup_info
  ) %>%
  addLegend(
    position = "bottomright",
    title = "Instalações Eólicas",
    colors = "darkred",
    labels = "Usinas"
  )

2. Teste 2 – Bird movement and tide height

Introdução

Este relatório apresenta a análise de movimentos de aves em relação à altura da maré, com foco em responder o Teste 2 de avaliação de dados e visualização temporal e espacial.
O objetivo é explorar padrões de comportamento de aves, sua distribuição altitudinal e correlação com fatores ambientais, como o ciclo das marés.

Perguntas do Teste 2

As análises abaixo foram desenvolvidas para responder às seguintes questões:

  1. Descrever o número de movimentos de aves por hora.
  2. Descrever a altura média de voo por hora.
  3. Verificar diferenças no número médio de movimentos de aves por hora considerando períodos diurnos e noturnos.
  4. Verificar diferenças na altura média de voo por hora considerando períodos diurnos e noturnos.
  5. Descrever o padrão de movimento das aves em relação à altura da maré e à distância para a próxima maré alta.
    • Sugestão: considerar a primeira observação de cada movimento para calcular a distância até a próxima maré alta.

2.1. Carregamento dos Dados e Pacotes

# ===============================
# Setup e carregamento de dados
# ===============================
library(tidyverse)
library(lubridate)
library(sf)
library(ggplot2)

# Carregar os dados do Teste 2
load("data/bird_track.rda")   # dataset de movimentos de aves
load("data/tidepeak.rda")      # dataset de altura da maré

# Visualização inicial
head(bird_track)
##                 date_l                             track_id speed heading
##                 <POSc>                               <char> <int>   <int>
## 1: 2021-03-13 21:51:59 4463cc43-a679-46d4-af30-b6886cf8c8cc     2      89
## 2: 2021-03-13 21:52:02 4463cc43-a679-46d4-af30-b6886cf8c8cc     2      89
## 3: 2021-03-13 21:52:02 4b4d3c6f-ee2a-42c0-9d37-6f427d00b72d    12     125
## 4: 2021-03-13 21:52:04 4b4d3c6f-ee2a-42c0-9d37-6f427d00b72d    13     122
## 5: 2021-03-13 21:52:04 818f11ea-c7b5-4f78-b434-752b92702930     2      84
## 6: 2021-03-13 21:52:04 9a42b4f5-8338-4d11-a614-5eed65c48349    15     247
##        altm
##       <num>
## 1:  41.7576
## 2:  41.7576
## 3: 296.5704
## 4: 288.0360
## 5:  81.9912
## 6:  81.9912
head(tidepeak)
##                   data      mare altura             dateref
##                 <POSc>    <char>  <num>              <POSc>
## 1: 2020-01-01 00:27:00 Baixa-Mar   1.50 2020-01-01 00:27:00
## 2: 2020-01-01 07:00:00 Preia-Mar   3.44 2020-01-01 07:00:00
## 3: 2020-01-01 13:05:00 Baixa-Mar   1.45 2020-01-01 13:05:00
## 4: 2020-01-01 19:32:00 Preia-Mar   3.15 2020-01-01 19:32:00
## 5: 2020-01-02 01:15:00 Baixa-Mar   1.64 2020-01-02 01:15:00
## 6: 2020-01-02 07:46:00 Preia-Mar   3.30 2020-01-02 07:46:00
# Definir opções globais do chunk
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

2.2 Limpeza, ETL e Padronização dos Dados

Nesta etapa foi realizada a preparação dos dados de monitoramento de aves (bird_track) e dos dados de maré (tidepeak). O processo incluiu:

  1. Remoção de valores ausentes;

  2. Filtragem de altitudes inválidas;

  3. Padronização de formatos de data e hora;

  4. Definição consistente de fuso horário (Europe/Lisbon);

  5. Ordenação dos dados para análises temporais;

  6. Criação de variáveis auxiliares (hora, data e período dia/noite).

OBS.: A padronização temporal é especialmente importante para garantir a correta integração entre movimentos de aves e os dados de maré.

# ==========================================
# 1. Exploração inicial da estrutura
# ==========================================

glimpse(bird_track)
## Rows: 79,453
## Columns: 5
## $ date_l   <dttm> 2021-03-13 21:51:59, 2021-03-13 21:52:02, 2021-03-13 21:52:0…
## $ track_id <chr> "4463cc43-a679-46d4-af30-b6886cf8c8cc", "4463cc43-a679-46d4-a…
## $ speed    <int> 2, 2, 12, 13, 2, 15, 1, 59, 12, 17, 2, 14, 1, 59, 12, 18, 4, …
## $ heading  <int> 89, 89, 125, 122, 84, 247, 90, 96, 122, 92, 84, 251, 90, 96, …
## $ altm     <dbl> 41.7576, 41.7576, 296.5704, 288.0360, 81.9912, 81.9912, 41.45…
glimpse(tidepeak)
## Rows: 2,825
## Columns: 4
## $ data    <dttm> 2020-01-01 00:27:00, 2020-01-01 07:00:00, 2020-01-01 13:05:00…
## $ mare    <chr> "Baixa-Mar", "Preia-Mar", "Baixa-Mar", "Preia-Mar", "Baixa-Mar…
## $ altura  <dbl> 1.50, 3.44, 1.45, 3.15, 1.64, 3.30, 1.57, 3.06, 1.75, 3.18, 1.…
## $ dateref <dttm> 2020-01-01 00:27:00, 2020-01-01 07:00:00, 2020-01-01 13:05:00…
# ==========================================
# 2. Remoção de valores ausentes e inconsistentes
# ==========================================

n_before <- nrow(bird_track)

bird_track_clean <- bird_track %>%
  filter(
    !is.na(track_id),
    !is.na(date_l),
    !is.na(altm),
    altm >= 0
  )

n_after <- nrow(bird_track_clean)

cat("Rows removed:", n_before - n_after, "\n")
## Rows removed: 336
# Limpeza dados de maré
tidepeak_clean <- tidepeak %>%
  filter(
    !is.na(data),
    !is.na(altura)
  )

# ==========================================
# 3. Padronização de classes temporais
# ==========================================

str(bird_track_clean$date_l)
##  POSIXct[1:79117], format: "2021-03-13 21:51:59" "2021-03-13 21:52:02" "2021-03-13 21:52:02" ...
str(tidepeak$data)
##  POSIXct[1:2825], format: "2020-01-01 00:27:00" "2020-01-01 07:00:00" "2020-01-01 13:05:00" ...
str(tidepeak$dateref)
##  POSIXct[1:2825], format: "2020-01-01 00:27:00" "2020-01-01 07:00:00" "2020-01-01 13:05:00" ...
# Definir timezone oficial
tz_lisbon <- "Europe/Lisbon"

bird_track_clean$date_l <- as.POSIXct(bird_track_clean$date_l, tz = tz_lisbon)
tidepeak_clean$data <- as.POSIXct(tidepeak_clean$data, tz = tz_lisbon)
tidepeak_clean$dateref <- as.POSIXct(tidepeak_clean$dateref, tz = tz_lisbon)

# Conferir timezone
attr(bird_track_clean$date_l, "tzone")
## [1] "Europe/Lisbon"
attr(tidepeak_clean$data, "tzone")
## [1] "Europe/Lisbon"
# ==========================================
# 4. Ordenação temporal
# ==========================================

bird_track_clean <- bird_track_clean %>%
  arrange(track_id, date_l)

tidepeak_clean <- tidepeak_clean %>%
  arrange(data)

# ==========================================
# 5. Criação de variáveis auxiliares
# ==========================================

bird_track_clean <- bird_track_clean %>%
  mutate(
    hour = lubridate::hour(date_l),
    date = as.Date(date_l),
    period = ifelse(hour >= 6 & hour < 18, "Day", "Night")
  )

2.3 Descreva o número de movimentos de aves por hora.

Observa-se uma variação expressiva no número de movimentos ao longo das horas do dia. Durante o período noturno, os maiores valores concentram-se entre 19h e 20h, com destaque para 19h, que apresenta o pico máximo de 28.429 movimentos. Já nas primeiras horas da madrugada (0h–5h), os valores permanecem elevados, especialmente às 5h (6.847 movimentos).

No período diurno, os registros são mais heterogêneos e, em geral, apresentam valores inferiores aos observados à noite, com exceção das 9h (3.102 movimentos). Em alguns horários, como 11h, o número de movimentos é bastante reduzido (apenas 6 registros).

Esses resultados indicam uma predominância de atividade no período noturno, com picos marcantes no início da noite. O expressivo aumento observado às 19h pode estar associado ao ritmo circadiano das aves, particularmente ao comportamento de deslocamento pré-empoleiramento (pre-roosting movements), quando os indivíduos intensificam seus movimentos em direção aos dormitórios. Esse padrão é frequentemente descrito como parte da organização temporal da atividade diária, refletindo a transição do período de forrageamento para o estabelecimento em áreas de abrigo coletivo ou individual, em resposta à redução da luminosidade e a fatores de segurança antipredatória.

# Número total de movimentos por hora e período
mov_hour <- bird_track_clean[
  , .(n_movements = .N),
  by = .(hour, period)
][order(hour, period)]

mov_hour
##      hour period n_movements
##     <int> <char>       <int>
##  1:     0  Night        3637
##  2:     1  Night        2830
##  3:     2  Night        2565
##  4:     3  Night        3751
##  5:     4  Night        5554
##  6:     5  Night        6847
##  7:     9    Day        3102
##  8:    10    Day         200
##  9:    11    Day           6
## 10:    12    Day        1272
## 11:    13    Day         345
## 12:    17    Day         487
## 13:    18  Night        1365
## 14:    19  Night       28429
## 15:    20  Night       12574
## 16:    21  Night         397
## 17:    22  Night        3284
## 18:    23  Night        2472
# Plot
ggplot(mov_hour, aes(x = hour, y = n_movements)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = 0:23) +
  labs(
    title = "Number of Bird Movements per Hour",
    x = "Hour of the Day",
    y = "Number of Movements"
  ) +
  theme_minimal()

2.4 Altura média de voo por hora

Observa-se uma variação clara na altura média de voo ao longo do dia. Durante a madrugada (0h–5h), as aves mantêm alturas relativamente estáveis, variando aproximadamente entre 129 m e 154 m, com leve elevação às 2h (153,66 m) e às 5h (149,63 m).

No período diurno, ocorre um aumento expressivo na altura média de voo, especialmente entre 11h e 13h, quando são registrados os maiores valores do dia, com destaque para 11h (282,96 m), representando o pico máximo observado. Esse padrão sugere possível influência de fatores atmosféricos, como a formação de térmicas e maior instabilidade convectiva ao longo do final da manhã e início da tarde, que favorecem voos em maiores altitudes.

No início da noite (18h–20h), as alturas permanecem elevadas (entre 215 m e 243 m), podendo refletir deslocamentos direcionais antes do período de repouso. Após as 21h, observa-se novamente uma redução gradual da altura média, retornando a valores próximos aos registrados durante a madrugada.

De forma geral, o padrão indica maior utilização da coluna vertical do espaço aéreo durante o período central do dia, possivelmente associada a condições ambientais favoráveis ao ganho de altitude, enquanto durante a madrugada predominam voos em altitudes mais moderadas e estáveis.

# Altura média de voo por hora.
height_hour <- bird_track_clean[
  , .(mean_height = mean(altm, na.rm = TRUE)),
  by = hour
][order(hour)]

height_hour
##      hour mean_height
##     <int>       <num>
##  1:     0    130.5246
##  2:     1    137.4580
##  3:     2    153.6632
##  4:     3    129.3559
##  5:     4    130.9728
##  6:     5    149.6286
##  7:     9    115.9674
##  8:    10    126.4707
##  9:    11    282.9560
## 10:    12    249.2936
## 11:    13    246.1547
## 12:    17    179.7131
## 13:    18    215.8319
## 14:    19    220.7991
## 15:    20    243.0086
## 16:    21    183.7053
## 17:    22    154.6670
## 18:    23    157.2254
# Plot
ggplot(height_hour, aes(x = hour, y = mean_height)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = 0:23) +
  labs(
    title = "Mean Flight Height per Hour",
    x = "Hour of Day",
    y = "Mean Flight Height (m)"
  ) +
  theme_minimal()

2.5 Diferenças no número médio de movimentos por hora entre dia e noite

A comparação entre os períodos diurno e noturno revelou diferenças marcantes na intensidade de movimentos por hora.

No período noturno, a média foi de 6.142 movimentos/hora (DP = 7.699; n = 12 horas), enquanto no período diurno a média foi substancialmente inferior, com 902 movimentos/hora (DP = 1.162; n = 6 horas). Essa diferença descritiva já sugere uma predominância expressiva de atividade durante a noite.

A normalidade dos dados foi avaliada por meio do teste de Shapiro–Wilk, adotando-se um nível de significância de 5% (IC de 95%), que corresponde ao padrão dos testes estatísticos utilizados. Os resultados indicaram violação do pressuposto de normalidade em ambos os períodos:

Dia: W = 0,786; p = 0,0437

Noite: W = 0,655; p = 0,00032

Como os valores de p foram inferiores a 0,05, rejeita-se a hipótese nula de normalidade ao nível de 95% de confiança.

Diante dessa assimetria, aplicou-se o teste não paramétrico de Wilcoxon (rank sum test), também considerando α = 0,05 (IC 95%), para comparar os períodos. O teste indicou diferença estatisticamente significativa entre dia e noite (W = 7; p = 0,00474), confirmando que o número de movimentos por hora é significativamente maior durante o período noturno.

Ecologicamente, esse padrão reforça a evidência de maior atividade noturna, possivelmente associada a deslocamentos crepusculares e movimentos direcionados a áreas de dormitório, o que pode explicar os picos observados no início da noite.

# Garantir que period é fator
mov_hour[, period := factor(period, levels = c("Day", "Night"))]

# Estatísticas descritivas
mean_mov_period <- mov_hour[
  , .(
    mean_movements = mean(n_movements),
    sd_movements   = sd(n_movements),
    n_hours        = .N
  ),
  by = period
]

mean_mov_period
##    period mean_movements sd_movements n_hours
##    <fctr>          <num>        <num>   <int>
## 1:  Night       6142.083     7699.715      12
## 2:    Day        902.000     1162.325       6
# Visualização: violino + boxplot
ggplot(mov_hour, aes(x = period, y = n_movements, fill = period)) +
  geom_violin(trim = FALSE, alpha = 0.4, color = NA) +
  geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.08, alpha = 0.5, size = 2) +
  scale_fill_manual(values = c("Day" = "#FDB813",
                               "Night" = "#1F3B73")) +
  labs(
    title = "Movimentos de Aves por Hora: Dia vs Noite",
    x = "Período",
    y = "Número de Movimentos por Hora"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Histograma
ggplot(mov_hour, aes(x = n_movements)) +
  geom_histogram(bins = 10, fill = "grey70", color = "white") +
  facet_wrap(~ period, scales = "free") +
  theme_minimal() +
  labs(title = "Distribuição dos Movimentos por Hora")

# QQ-plot
ggplot(mov_hour, aes(sample = n_movements)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~ period)

# Teste de normalidade
shapiro.test(mov_hour$n_movements[mov_hour$period == "Day"])
## 
##  Shapiro-Wilk normality test
## 
## data:  mov_hour$n_movements[mov_hour$period == "Day"]
## W = 0.78592, p-value = 0.04374
shapiro.test(mov_hour$n_movements[mov_hour$period == "Night"])
## 
##  Shapiro-Wilk normality test
## 
## data:  mov_hour$n_movements[mov_hour$period == "Night"]
## W = 0.65517, p-value = 0.0003223
# Teste não paramétrico
wilcox.test(n_movements ~ period, data = mov_hour)
## 
##  Wilcoxon rank sum exact test
## 
## data:  n_movements by period
## W = 7, p-value = 0.00474
## alternative hypothesis: true location shift is not equal to 0

2.6 Diferenças na altitude média de voo por hora entre dia e noite

A análise da altitude média de voo por hora indica que, descritivamente, as aves tendem a voar em altitudes mais elevadas durante o período diurno. A média das altitudes horárias durante o dia foi de 200,09 m (DP = 69,74; n = 6 horas), enquanto no período noturno foi de 167,24 m (DP = 39,29; n = 12 horas). Além da diferença nas médias, observa-se maior variabilidade das altitudes durante o dia.

A normalidade dos dados foi avaliada pelo teste de Shapiro–Wilk, adotando nível de significância de 5% (IC de 95%). Para o período diurno, não houve evidência de violação da normalidade (W = 0,897; p = 0,357). Entretanto, para o período noturno, a normalidade foi rejeitada (W = 0,851; p = 0,038), indicando distribuição não normal.

Diante da assimetria observada em pelo menos um dos grupos, foi aplicado o teste não paramétrico de Wilcoxon (rank sum test) para comparar as altitudes médias horárias entre os períodos. O resultado não indicou diferença estatisticamente significativa entre dia e noite (W = 44; p = 0,494; IC de 95%).

Assim, embora haja uma diferença descritiva nas médias — com valores mais elevados durante o dia — essa variação não é estatisticamente significativa ao nível de 5%, sugerindo que a altitude média de voo por hora não difere de forma consistente entre os períodos diurno e noturno dentro do conjunto de dados analisado.

# Calcular altura média por hora com período
height_hour <- bird_track_clean[
  , .(mean_height = mean(altm, na.rm = TRUE)),
  by = .(hour, period)
][order(hour)]

height_hour
##      hour period mean_height
##     <int> <char>       <num>
##  1:     0  Night    130.5246
##  2:     1  Night    137.4580
##  3:     2  Night    153.6632
##  4:     3  Night    129.3559
##  5:     4  Night    130.9728
##  6:     5  Night    149.6286
##  7:     9    Day    115.9674
##  8:    10    Day    126.4707
##  9:    11    Day    282.9560
## 10:    12    Day    249.2936
## 11:    13    Day    246.1547
## 12:    17    Day    179.7131
## 13:    18  Night    215.8319
## 14:    19  Night    220.7991
## 15:    20  Night    243.0086
## 16:    21  Night    183.7053
## 17:    22  Night    154.6670
## 18:    23  Night    157.2254
# Estatísticas descritivas por período
mean_height_period <- height_hour[
  , .(
    mean_height = mean(mean_height),
    sd_height   = sd(mean_height),
    n_hours     = .N
  ),
  by = period
]

mean_height_period
##    period mean_height sd_height n_hours
##    <char>       <num>     <num>   <int>
## 1:  Night    167.2367  39.28925      12
## 2:    Day    200.0926  69.74434       6
# Visualização (violino + boxplot)
ggplot(height_hour, aes(x = period, y = mean_height, fill = period)) +
  geom_violin(trim = FALSE, alpha = 0.4, color = NA) +
  geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.08, alpha = 0.5, size = 2) +
  scale_fill_manual(values = c("Day" = "#FDB813",
                               "Night" = "#1F3B73")) +
  labs(
    title = "Mean Flight Height per Hour: Day vs Night",
    x = "Period",
    y = "Mean Flight Height (m)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Verificar normalidade (IC 95%)
shapiro.test(height_hour[period == "Day"]$mean_height)
## 
##  Shapiro-Wilk normality test
## 
## data:  height_hour[period == "Day"]$mean_height
## W = 0.89709, p-value = 0.357
shapiro.test(height_hour[period == "Night"]$mean_height)
## 
##  Shapiro-Wilk normality test
## 
## data:  height_hour[period == "Night"]$mean_height
## W = 0.85121, p-value = 0.03799
# QQ-plot
ggplot(height_hour, aes(sample = mean_height)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~ period)

# Teste estatístico não paramétrico (IC 95%)
wilcox.test(mean_height ~ period, data = height_hour)
## 
##  Wilcoxon rank sum exact test
## 
## data:  mean_height by period
## W = 44, p-value = 0.4936
## alternative hypothesis: true location shift is not equal to 0

Desafio 2 - Descrever o padrão de movimento das aves em relação à altura da maré e à distância até a próxima preia-mar

Nesta etapa, o objetivo é investigar como os movimentos das aves se relacionam com a altura da maré e com o tempo restante até a próxima preia-mar, considerando apenas o primeiro registro de cada movimento (track), conforme sugerido no enunciado.

Descrição detalhada do procedimento adotado

Inicialmente, os dados de rastreamento (bird_track_clean) são ordenados por track_id e date_l, garantindo que as observações estejam organizadas cronologicamente dentro de cada movimento. Em seguida, é utilizada a operação:

first_mov <- bird_track_clean[
  , .SD[1],
  by = track_id
]

Aqui, .SD[1] seleciona a primeira observação de cada track_id, ou seja, o ponto inicial de cada movimento. Essa etapa é fundamental para evitar pseudorreplicação, já que múltiplos registros pertencentes ao mesmo movimento poderiam enviesar a análise temporal em relação à maré.

Posteriormente, os dados de maré (tidepeak) são filtrados para manter apenas eventos de preia-mar:

tide_high <- tidepeak[mare == "Preia-Mar"]

Os dados são então ordenados cronologicamente pela variável dateref, que representa o horário de referência da maré alta.

Para relacionar cada movimento com a maré, é realizado um join temporal não exato (rolling join):

first_mov_tide <- tide_high[
  first_mov,
  on = .(dateref = date_l),
  roll = Inf
]

Essa operação associa cada primeiro registro de movimento ao evento de preia-mar mais próximo no tempo (de acordo com a ordenação), permitindo vincular o movimento ao ciclo de maré correspondente.

Em seguida, é calculada a variável:

first_mov_tide[, hours_to_high_tide :=
                 as.numeric(difftime(dateref, data, units = "hours"))
]

Essa variável representa a distância temporal (em horas) até a próxima preia-mar, sendo um indicador contínuo da posição do movimento dentro do ciclo de maré.

Relação entre movimentos e altura da maré

A distribuição dos movimentos em função da altura da maré revela uma forte concentração em determinados níveis. Observa-se que:

A altura 3.55 m concentra aproximadamente 46,9% dos movimentos.

A altura 3.8 m representa cerca de 29,0%.

A altura 3.99 m corresponde a 17,6%.

Alturas menores (3.4 m e 3.51 m) apresentam proporções significativamente inferiores.

# Transformar variáveis em fatores para análise categórica
first_mov_tide$altura_f <- as.factor(first_mov_tide$altura)
first_mov_tide$hours_int <- as.factor(floor(first_mov_tide$hours_to_high_tide))

# ----------------------------
# Movimento vs Altura da Maré
# ----------------------------

ggplot(first_mov_tide, aes(x = altura_f)) +
  geom_bar() +
  theme_minimal() +
  labs(
    title = "Movimentos em função da altura da maré",
    x = "Altura da Maré (m)",
    y = "Frequência de Movimentos"
  )

# Tabela de contagem absoluta
table(first_mov_tide$altura_f)
## 
##  3.4 3.51 3.55  3.8 3.99 
##   54 1087 8212 5073 3076
# Frequência relativa (%)
prop.table(table(first_mov_tide$altura_f)) * 100
## 
##        3.4       3.51       3.55        3.8       3.99 
##  0.3085362  6.2107188 46.9203520 28.9852588 17.5751343

Relação entre movimentos e distância até a próxima preia-mar

Ao analisar a distância temporal até a próxima preia-mar, observa-se:

O maior percentual de movimentos ocorre 1 hora antes da preia-mar (32,4%).

Em seguida, destaca-se o intervalo de 2 horas antes (22,4%).

Movimentos exatamente no momento da maré alta (0 h) representam 7,5%.

Após 3–4 horas, há queda acentuada na frequência.

A mediana da distância temporal é 2,24 horas, enquanto a média é 4,14 horas, indicando leve assimetria à direita (presença de valores mais distantes da maré alta).

Esses resultados sugerem que os movimentos tendem a se concentrar nas horas imediatamente anteriores à preia-mar, possivelmente refletindo sincronização comportamental com o ciclo de maré.

# -------------------------------------------
# Movimento vs Distância até Preia-Mar (horas)
# -------------------------------------------

ggplot(first_mov_tide, aes(x = hours_int)) +
  geom_bar() +
  theme_minimal() +
  labs(
    title = "Movimentos em função da distância até a próxima preia-mar",
    x = "Horas até a próxima preia-mar",
    y = "Frequência de Movimentos"
  )

# Contagem absoluta
table(first_mov_tide$hours_int)
## 
##    0    1    2    3    4    6    7    8    9   10   11   12 
## 1306 5672 3925 1059  163  812  888  909  936  840  857  135
# Frequência relativa (%)
prop.table(table(first_mov_tide$hours_int)) * 100
## 
##          0          1          2          3          4          6          7 
##  7.4620043 32.4077248 22.4260085  6.0507371  0.9313221  4.6394698  5.0737059 
##          8          9         10         11         12 
##  5.1936921  5.3479602  4.7994515  4.8965832  0.7713404
# Estatísticas descritivas
summary(first_mov_tide$hours_to_high_tide)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.002222  1.690278  2.236250  4.142357  7.067847 12.231944