In pratikunterwegs/holeybirds: Link bird moult with movement.

Vegetation, visibility, and land use in the study area

Load libraries

# libs for data
library(data.table)
library(glue)
library(dplyr)

library(sf)
library(terra)
# to plot rasters
library(stars)

library(ggplot2)
library(colorspace)
library(patchwork)

Load data

# load rasters using terra
ndvi <- rast("data/rasters/raster_hula_ndvi_2039.tif")
vis <- rast("data/rasters/raster_hula_visibility.tif")

# load land use as sf
landcover <- st_read("data/spatial/hula_lc_vector")
landcover <- st_crop(landcover, vis)

# crop ndvi by vis
ndvi <- crop(ndvi, vis)

NDVI across the landscape

ndvi <- st_as_stars(ndvi)

fig_ndvi <-
  ggplot() +
  geom_stars(
    data = ndvi,
    # downsample = 5
  ) +
  coord_sf(
    crs = 2039,
    expand = F
  ) +
  scale_x_continuous(
    breaks = seq(35.58, 35.64, length.out = 3)
  ) +
  ggspatial::annotation_scale() +
  scico::scale_fill_scico(
    direction = -1,
    palette = "bamako",
    limits = c(0.05, 0.7),
    breaks = c(0.05, seq(0.3, 0.7, 0.2)),
    name = "NDVI",
    na.value = "aliceblue"
  ) +
  theme_test() +
  theme(
    legend.position = "bottom",
    legend.key.height = unit(2, "mm"),
    legend.title = element_text(
      vjust = 1
    ),
    axis.title = element_blank(),
    axis.text = element_text(
      size = 8
    ),
    axis.text.y = element_text(
      angle = 90, hjust = 0.5
    )
  )

Long term NDVI

# load files from csv
ndvi_mean_longterm <- fread("data/rasters/ndvi_long_term/ndvi_week_mean_20100101.csv")
ndvi_sd_longterm <- fread("data/rasters/ndvi_long_term/ndvi_week_sd_20100101.csv")

ndvi_stats <- Map(
  list(ndvi_mean_longterm, ndvi_sd_longterm), list("ndvi_mean", "ndvi_sd"),
  f = function(df, val) {
    melt(
      df[, -c("system:index")],
      id = c("OBJECTID", "SHAPE_Leng", "SHAPE_Area", "Name"),
      variable = "date", value = val
    )[, date := as.Date(stringr::str_extract(date, "\\d+"), format = "%Y%m%d")]
  }
)
ndvi_stats <- Reduce(x = ndvi_stats, f = merge.data.table)

# save NDVI stats and track them
fwrite(ndvi_stats, "data/results/ndvi_stats_longterm.csv")

# summary ndvi stats by landcover and weighted by area
# remove water areas
ndvi_stats_summary <- ndvi_stats[, week := week(date)][, list(
  ndvi_mean = weighted.mean(ndvi_mean, SHAPE_Area, na.rm = TRUE),
  ndvi_sd = weighted.mean(ndvi_sd, SHAPE_Area, na.rm = TRUE),
  # artificial date for plotting
  date = as.Date(paste(2016, week, 1, sep = "-"), "%Y-%U-%u")
), by = c("Name", "week")]
ndvi_stats_summary <- ndvi_stats_summary[!Name %in% c("W", "C"), ]

# ndvi stats for 2016
ndvi_stats_2016 <- ndvi_stats[year(date) == 2016, ]
ndvi_stats_2016 <- ndvi_stats_2016[, week := week(date)][, list(
  ndvi_mean = weighted.mean(ndvi_mean, SHAPE_Area, na.rm = TRUE),
  ndvi_sd = weighted.mean(ndvi_sd, SHAPE_Area, na.rm = TRUE),
  # artificial date for plotting
  date = as.Date(paste(2016, week, 1, sep = "-"), "%Y-%U-%u")
), by = c("Name", "week")]
ndvi_stats_2016 <- ndvi_stats_2016[!Name %in% c("W", "C"), ]

fwrite(ndvi_stats_2016, "data/results/ndvi_stats_2016_summary.csv")

# Read in study duration
id_data <- fread("data/results/data_tracking_metrics_rrv.csv")
min_tracking <- min(id_data$start_tracking)
max_tracking <- max(id_data$end_tracking)

id_data <- id_data[, list(
  min_tracking = min(start_tracking),
  max_tracking = max(start_tracking)
),
by = c("sp", "treat")
]

Plot for supplementary material.

id_data <- split(id_data[sp != "Hirundo", ], by = "sp")

plots_ndvi_tracking <- Map(
  id_data, names(id_data),
  f = function(df, sp) {
    ggplot() +
      geom_ribbon(
        data = ndvi_stats_summary,
        aes(
          date,
          ymin = ndvi_mean - ndvi_sd,
          ymax = ndvi_mean + ndvi_sd,
          group = Name
        ),
        fill = "grey90"
      ) +
      geom_line(
        data = ndvi_stats_summary,
        aes(
          date, ndvi_mean,
          group = Name
        ),
        colour = "grey"
      ) +
      geom_pointrange(
        data = ndvi_stats_2016,
        aes(
          x = date, y = ndvi_mean,
          ymin = ndvi_mean - ndvi_sd,
          ymax = ndvi_mean + ndvi_sd,
          group = Name
        ),
        colour = "grey60",
        stroke = 0.2
      ) +
      geom_rect(
        data = df,
        aes(
          xmin = as.Date(min_tracking),
          xmax = as.Date(max_tracking),
          ymin = -1, ymax = 1,
          fill = treat,
          col = treat
        ),
        alpha = 0.6
      ) +
      scale_x_date(
        date_breaks = "2 month",
        date_labels = "%b"
      ) +
      scale_fill_discrete_sequential(
        name = NULL,
        palette = "Batlow",
        l1 = 30, l2 = 60,
        breaks = c("NonMoulting", "Moulting", "Manipulated"),
        labels = c("Non-molting", "Molting", "Manipulated")
      ) +
      scale_colour_discrete_sequential(
        name = NULL,
        palette = "Batlow",
        l1 = 50, l2 = 50,
        guide = "none"
      ) +
      facet_grid(
        rows = vars(treat),
        cols = vars(Name),
        labeller = labeller(
          Name = c(
            "B" = "Settlements",
            "O" = "Open areas",
            "R" = "Reedbeds",
            "T" = "Trees"
          )
        )
      ) +
      coord_cartesian(
        xlim = as.Date(c("2016-06-01", "2016-12-01")),
        ylim = c(0, 0.6),
        expand = TRUE
      ) +
      labs(
        x = NULL,
        y = "NDVI",
        title = sp
      ) +
      theme_test(
        base_size = 10,
        base_family = "Arial"
      ) +
      theme(
        legend.position = "top",
        strip.background = element_blank(),
        title = element_text(
          face = "italic"
        )
      )
  }
)

Map(
  plots_ndvi_tracking, names(plots_ndvi_tracking),
  f = function(gg, sp) {
    ggsave(gg,
      filename = sprintf("figures/fig_tracking_ndvi_%s.png", sp),
      height = 5, width = 6
    )
  }
)

Visibility across the landscape

vis <- st_as_stars(vis)

fig_vis <-
  ggplot() +
  geom_stars(
    data = vis,
    downsample = 3
  ) +
  coord_sf(
    crs = 2039,
    expand = F
  ) +
  scale_x_continuous(
    breaks = seq(35.58, 35.64, length.out = 3)
  ) +
  ggspatial::annotation_scale() +
  scico::scale_fill_scico(
    direction = 1,
    palette = "davos",
    limits = c(0, 1),
    breaks = seq(0., 1, 0.2),
    name = "Vis. index"
  ) +
  theme_test() +
  theme(
    legend.position = "bottom",
    legend.key.height = unit(2, "mm"),
    legend.title = element_text(
      vjust = 1
    ),
    axis.title = element_blank(),
    axis.text = element_text(
      size = 8
    ),
    axis.text.y = element_text(
      angle = 90, hjust = 0.5
    )
  )

Plot landcover

# handle water classes
landcover <- mutate(
  landcover,
  class = if_else(
    Name %in% c("C", "W"),
    "W", Name
  )
)

fig_landcover <-
  ggplot(landcover) +
  geom_sf(
    aes(
      fill = class
    ),
    col = "transparent"
  ) +
  scico::scale_fill_scico_d(
    palette = "hawaii",
    direction = 1,
    labels = c(
      "B" = "Settlements",
      "O" = "Open areas",
      "R" = "Reedbeds",
      "T" = "Trees",
      "W" = "Water"
    ),
    name = NULL
  ) +
  coord_sf(
    crs = 2039,
    expand = F
  ) +
  scale_x_continuous(
    breaks = seq(35.58, 35.64, length.out = 3)
  ) +
  ggspatial::annotation_scale() +
  theme_test() +
  theme(
    legend.position = "bottom",
    legend.key.height = unit(2, "mm"),
    legend.key.width = unit(2, "mm"),
    legend.title = element_text(
      vjust = 1
    ),
    legend.text = element_text(
      size = 6
    ),
    panel.background = element_rect(
      fill = "grey99"
    ),
    plot.background = element_blank(),
    axis.title = element_blank(),
    axis.text = element_text(
      size = 8
    ),
    axis.text.y = element_text(
      angle = 90, hjust = 0.5
    )
  ) +
  guides(
    fill = guide_legend(
      nrow = 2
    )
  )

Merge figures

fig_landscape_maps <- wrap_plots(
  fig_landcover, fig_ndvi, fig_vis
) + plot_annotation(
  tag_levels = "A"
) &
  theme(
    plot.tag = element_text(
      size = 10, face = "bold"
    )
  )

# save overall figure
ggsave(
  fig_landscape_maps,
  filename = "figures/fig_landscape_maps.png",
  height = 125, width = 200, units = "mm"
)

Model correlation NDVI and visibility

Reload rasters as terra objects.

# load rasters using terra
ndvi <- rast("data/rasters/raster_hula_ndvi_2039.tif")
vis <- rast("data/rasters/raster_hula_visibility.tif")

# make samples
samples <- st_sample(
  landcover,
  size = 10000,
  type = "hexagonal"
)

# extract data at samples
vis_samp <- terra::extract(
  vis,
  vect(samples)
)
setDT(vis_samp)

ndvi_samp <- extract(
  ndvi,
  vect(samples)
)
setDT(ndvi_samp)

# lc samples
lc_samp <- extract(
  vect(
    dplyr::select(landcover, Name)
  ),
  vect(samples)
)
setDT(lc_samp)

# merge data
env_data <- merge(vis_samp, ndvi_samp)
env_data <- merge(env_data, lc_samp, by.x = "ID", by.y = "id.x")

# filter for valid values
env_data <- env_data[!is.nan(raster_hula_visibility)]
env_data <- setnames(env_data[, !("ID")], c("vis", "ndvi", "lc"))

# landcover classifications
# B = settlement
# C = canal
# O = agriculture/open
# R = reedbeds
# T = trees
# W = water

# remove water
env_data <- env_data[!lc %in% c("C", "W")]

# set landcover factor levels
env_data$lc <- factor(
  env_data$lc,
  levels = c("T", "R", "B", "O")
)

Fit a simple GAM

# load library
library(mgcv)

# fit a bam for large data
mod <- bam(
  vis ~ s(ndvi, by = lc) +
    s(lc, bs = "re"),
  data = env_data
)

summary(mod)

model_summary <- summary(mod)

writeLines(
  capture.output(model_summary),
  con = "data/results/model_coefs_vis_ndvi.txt"
)

gratia::draw(mod, fixed = T)

Get predicted values

# make prediction table
pred_data <- CJ(
  lc = factor(unique(env_data$lc),
    levels = c("T", "R", "B", "O", "W")
  ),
  ndvi = seq(0, 1, 0.02)
)

# get prediction
pred <- predict(mod, newdata = pred_data, allow.new.levels = T, se.fit = T)

pred_data$pred <- pred$fit
pred_data$se <- pred$se.fit

fig_vis_ndvi <-
  ggplot(pred_data) +
  geom_bin_2d(
    data = env_data,
    aes(
      ndvi, vis,
      fill = ..count..
    ),
    show.legend = T
  ) +
  geom_ribbon(
    aes(
      ndvi,
      ymin = pred - se,
      ymax = pred + se,
      group = lc
    ),
    show.legend = F,
    size = 0.3,
    col = "grey",
    fill = "transparent"
  ) +
  geom_line(
    aes(
      ndvi, pred
    ),
    col = "indianred"
  ) +
  facet_wrap(
    facets = vars(lc),
    labeller = labeller(
      lc = c(
        "B" = "Settlements",
        "O" = "Open areas",
        "R" = "Reedbeds",
        "T" = "Trees"
      )
    )
  ) +
  labs(
    x = "NDVI",
    y = "Visibility index"
  ) +
  scico::scale_fill_scico(
    palette = "lapaz",
    direction = -1,
    name = "# Samples"
  ) +
  coord_cartesian(
    expand = T,
    ylim = c(0, 1),
    xlim = c(0, 0.8)
  ) +
  theme_test() +
  theme(
    legend.position = "bottom",
    # legend.justification = 1,
    legend.margin = margin(rep(0, 4)),
    legend.key.height = unit(2, "mm"),
    legend.key.width = unit(10, "mm"),
    strip.background = element_blank(),
    strip.text = element_text(
      face = "italic"
    )
  )

# save figure
ggsave(
  fig_vis_ndvi,
  filename = "figures/fig_spm_03.png",
  height = 90, width = 90, units = "mm"
)