Moon Cycle in relation to Sleeping site Selection

1 Libraries

We relied on a combination of spatial, statistical, and visualization packages in R to conduct this analysis. Each library played a specific role:

• tidyverse: A suite of packages for data manipulation and visualization. Core functions from dplyr, ggplot2, tibble, and readr were used to wrangle and plot data in a consistent and readable syntax.

• terra: Used for handling spatial raster data. We imported and processed the global canopy height raster (Lang et al. 2022), extracted values, and performed cropping and zooming for localized analysis.

• suncalc: Provided astronomical calculations, including moon illumination fraction, moonrise and moonset times, and phase information. This allowed us to incorporate lunar variables into our ecological questions.

• purrr: Part of the tidyverse ecosystem, purrr was essential for iterating over rows and applying functions (e.g., moon calculations) across date and location combinations in a clean and vectorized way.

• hexbin: Enabled the creation of hexagonal binned plots to explore the density distribution of sleeping sites in relation to canopy height and moonlight, which helps overcome overplotting in dense datasets.

• viridis: Provided perceptually uniform and colorblind-friendly color scales for raster and density plots, improving visual interpretation of ecological gradients.

• geosphere: Used for accurate geographic calculations such as generating random points within a radius from a known location (destPoint()), taking Earth’s curvature into account — crucial for ecological realism in spatial simulation.

library(tidyverse)
library(terra)
library(suncalc)
library(purrr)
library(hexbin)
library(viridis)
library(geosphere)
library(scales)

1.1 Description of the moon.rds dataset

The moon.rds file contains georeferenced data on sleeping sites of four groups of black capuchin monkeys (Sapajus nigritus) in **Iguazú National Park, Argentina, collected between 2011 and 2022.

Only night-time sleeping locations were considered. Specifically, each record corresponds to the last GPS point of the day (and always after sunset), collected after the monkeys had settled for the night, and does not include morning sleeping positions to avoid duplicates.

Each observation includes:

• GPS locations (longitude and latitude in WGS84), obtained either via:
• GPS collars, only including high-precision data based on HDOP and satellite number
• Direct field observations collected with handheld GPS devices
• The Date of the sleeping location record, which allows the estimation of lunar illumination (moon phases and altitude) for each night; detailed below.

This dataset is designed to support analyses exploring the relationship between lunar cycles and illumination and sleeping site selection, such as sleeping height or location in relation to moonlight.

moon <- readRDS("moon2.rds")

1.2 Calculating Lunar Illumination

To evaluate how lunar conditions might influence the sleeping site selection of black capuchin monkeys (Sapajus nigritus), we calculated key lunar variables using the suncalc package.

In the code below, for each observation date in the moon dataset, we extracted:

• moon_fraction: the proportion of the moon’s disc illuminated on that night (ranging from 0 for New Moon to 1 for Full Moon).
• moon_phase_angle: a normalized value (0 to 1) indicating the moon’s position within its cycle (e.g., 0 = New Moon, 0.5 = Full Moon).
• moon_altitude: which was setted for three specific hours during the night and then averaged, procedure is detailed below.

These lunar metrics might be ecologically relevant because moonlight intensity can influence nocturnal visibility and perceived predation risk. In primates, including capuchins, we predict that individuals may select higher sleeping sites during dark nights, which might be linked to a reduction in predation risk perception.

By incorporating lunar illumination data, we can test whether capuchin monkeys adjust their sleeping site height or location in response to varying moonlight conditions across nights and seasons.

moon <- moon %>%
  mutate(
    moon_fraction = sapply(fecha, function(f) getMoonIllumination(as.POSIXct(f))$fraction),
    moon_phase_angle = sapply(fecha, function(f) getMoonIllumination(as.POSIXct(f))$phase))

1.3 Classifying Moon Phase Categories

To further interpret lunar illumination patterns, we classified each observation into a moon phase category based on the moon_phase_angle value. The moon’s cycle is continuous, but for ecological analyses, it is often useful to group nights into broader phase categories that may reflect different levels of moonlight or perceived predation risk.

We defined the categories as follows:

• New Moon: moon_phase_angle < 0.1 or > 0.9
  (Very low illumination; darkest nights)
• Waxing: between 0.1 and 0.4
  (Increasing illumination; transition toward full moon)
• Full Moon: between 0.4 and 0.6
  (Peak illumination)
• Waning: between 0.6 and 0.9
  (Decreasing illumination after full moon)

Ecologically, these categories can help determine whether capuchin monkeys adjust their sleeping behavior, such as height in the canopy or site exposure, in response to the amount of ambient nocturnal light.

moon <- moon %>%
  mutate(
    moon_phase_name = case_when(
      moon_phase_angle < 0.1 | moon_phase_angle > 0.9 ~ "New Moon",
      moon_phase_angle >= 0.1 & moon_phase_angle < 0.4 ~ "Waxing",
      moon_phase_angle >= 0.4 & moon_phase_angle <= 0.6 ~ "Full Moon",
      moon_phase_angle > 0.6 & moon_phase_angle <= 0.9 ~ "Waning" ))

1.4 Moonrise and Moonset Times

To better understand the temporal availability of moonlight, we calculated the moonrise and moonset times for each observation using the getMoonTimes() function from the suncalc package. This helps determine whether the moon was actually present in the sky during the night — a crucial complement to the phase and illumination fraction.

We specified the location (latitude and longitude) of each observation, along with the observation date and the appropriate time zone for Iguazú National Park (America/Argentina/Cordoba)

These variables allow us to identify whether the moon was above the horizon during the active nighttime period. This is ecologically meaningful because a night with high lunar illumination on paper (e.g., a full moon) might still be functionally dark if the moon had already set by the time the monkeys were selecting or using their sleeping site.

In follow-up analyses, we can calculate:

Whether the moon was visible during key night hours (e.g., 19:00–05:00),

The duration the moon was above the horizon,

And how these metrics relate to sleeping site selection and height.

moon <- moon %>%
  mutate(
    moon_times = pmap(
      list(fecha, lat, lon),
      function(f, la, lo) {
        getMoonTimes(
          date = f,
          lat = la,
          lon = lo,
          keep = c("rise", "set"),
          tz = "America/Argentina/Cordoba"  
        )
      }
    )
  ) %>%
  mutate(
    moon_rise = map(moon_times, "rise") %>% flatten_dbl(),
    moon_set  = map(moon_times, "set") %>% flatten_dbl()
  ) %>%
  select(-moon_times)

moon <- moon %>%
  mutate(
    moon_rise = as.POSIXct(moon_rise, origin = "1970-01-01", tz = "America/Argentina/Cordoba"),
    moon_set  = as.POSIXct(moon_set,  origin = "1970-01-01", tz = "America/Argentina/Cordoba")
  )

2 Moon Visibility and Duration

We also calculated whether the moon was actually visible in the sky during each night, and for how long, based on the previously extracted moonrise and moonset times.

moon_duration_hr: the number of hours the moon was above the horizon that night.

moon_visible: a logical indicator (TRUE or FALSE) denoting whether the moon rose and set that night (i.e., was visible at some point).

These variables are ecologically important because they account for actual moonlight exposure at ground level. For example, even on a night with a full moon, if the moon had already set before the monkeys began selecting sleeping sites, its potential influence on behavior would be minimal. Conversely, on darker nights (low lunar fraction), a visible moon even for a few hours could still provide light that may influence perceived predation risk, vigilance, or site choice.

By integrating both moon phase and sky visibility, we can better capture the real nocturnal light environment experienced by the monkeys.

moon <- moon %>%
  mutate(
    moon_duration_hr = as.numeric(difftime(moon_set, moon_rise, units = "hours")),
    moon_visible = !is.na(moon_rise) & !is.na(moon_set) & moon_duration_hr > 0
  )

ggplot(moon, aes(x = group, y = moon_fraction, color = group)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.85, size = 2.8) +

  annotate("text", x = 4.5, y = 0, label = "🌑", size = 6) +
  annotate("text", x = 4.5, y = 0.5, label = "🌓", size = 6) +
  annotate("text", x = 4.5, y = 1, label = "🌕", size = 6) +

  scale_color_manual(values = c(
    "gun" = "darkorange3",
    "gue" = "midnightblue",
    "mac" = "goldenrod2",
    "spo" = "darkseagreen4"  #  "slategray4" could be a nice option
  )) +

  labs(
    title = "Observed Moon Illumination per Group",
    x = "Monkey Group",
    y = "Moon Illumination Fraction"
  ) +
  coord_cartesian(ylim = c(-0.05, 1.05)) +
  theme_bw(base_size = 14) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(face = "bold")
  )

Distribution of moon illumination fractions across monkey groups. Each point represents a recorded sleeping site, colored by group. Lunar phase icons on the right (🌑 new moon, 🌓 half moon, 🌕 full moon) provide a visual scale reference.

This plot shows the distribution of observed lunar illumination fractions (moon_fraction) across all recorded sleeping site events, separated by monkey group.

Each point represents a single observation, reflecting the proportion of the moon’s surface illuminated that night (0 = 🌑 New Moon, 1 = 🌕 Full Moon). The data appears to be well-distributed across the lunar cycle, with no strong evidence of sampling bias across groups — all groups show a broad range of moonlight conditions.

This preliminary exploration supports the idea that the dataset is suitable for investigating whether capuchin monkeys adjust their sleeping site choices (e.g., height or exposure) in response to nocturnal moonlight levels.

To explore the seasonality in sleeping site selection, two visualizations were created using the moon illumination fraction as the key variable. The first plot (circular) positions each observation according to the day of the year (angle) and moon illumination (radius), with a continuous color scale indicating the brightness of the moon. This layout helps reveal potential seasonal concentrations and lunar phase patterns in spatial behavior. The second plot, shown in a cartesian format, displays the same information across calendar months, allowing for clearer month-to-month comparisons. In both cases, the use of color to represent moonlight intensity provides an intuitive visual cue for interpreting temporal patterns in sleeping site usage.

moon <- moon %>%
  mutate(
    yday = yday(fecha),
    angle = 2 * pi * yday / 365 )

ggplot(moon, aes(x = angle, y = moon_fraction, fill = moon_fraction)) +
  geom_tile(width = 0.75, height = 0.19) +  # cuadraditos circulares

  scale_fill_viridis_c(
    option = "B",
    name = "Moon Illum.",
    guide = guide_colorbar(barheight = 8)) +

  coord_polar(start = -pi/2) +

  scale_x_continuous(
    breaks = seq(0, 2*pi, length.out = 13)[-13],
    labels = month.abb
  ) +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +

  labs(
    title = "Seasonality of Sleeping Sites by Moon Phase",
    x = NULL, y = NULL) +

  theme_minimal(base_size = 15) +
  theme(
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_blank(),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    legend.title = element_text(),
    plot.title = element_text(face = "bold")
  )

Circular visualization of sleeping site observations throughout the year. Each tile corresponds to a record, positioned by day of year (angle) and moon illumination (radius and color). This layout highlights seasonal patterns under varying moonlight conditions.

moon <- moon %>%
  mutate(month = month(fecha, label = TRUE, abbr = TRUE))

ggplot(moon, aes(x = month, y = moon_fraction, fill = moon_fraction)) +
  geom_violin(scale = "width", fill = "gray90", color = NA, alpha = 0.5) +
  geom_jitter(shape = 21, color = "black", size = 2.5, stroke = 0.3, width = 0.2, alpha = 0.8) +
  scale_fill_viridis_c(
  option = "B",
  name = "Moon Phase",
  limits = c(0, 1),  
  breaks = c(0, 0.25, 0.5, 0.75, 1),
 labels = c("0.00 - 🌑", "0.25 - 🌘", "0.50 - 🌗", "0.75 - 🌒", "1.00 - 🌕"),
  guide = guide_colorbar(barheight = 12)) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    title = "Seasonality of Sleeping Sites and Moon Illumination",
    x = "Month",
    y = "Moon Illumination Fraction"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(face = "bold"),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 13),
    legend.position = "right" )

Distribution of moonlight conditions during sleeping site observations throughout the year. Each point corresponds to an observation, jittered for visibility, and colored by lunar illumination. Violin plots represent the monthly density of observations. Emojis in the legend indicate approximate moon phases (new to full).

3 Canopy Height

We used the high-resolution global canopy height dataset produced by Lang et al. (2022), derived from GEDI LiDAR observations and machine learning models, to extract vegetation structure information across the Iguazú National Park. The original raster (ETH_GlobalCanopyHeight_10m_2020_S27W057) provides canopy height estimates at a 10-meter spatial resolution, allowing for detailed landscape-scale ecological analysis.

In this section, we cropped the raster to focus on the Iguazú National Park area (between -55 and -54 longitude, and -26.5 to -25.3 latitude), and applied a vegetation-like green color gradient to better visualize canopy height variations. This visualization enables us to interpret habitat vertical structure and its potential influence on monkey sleeping site selection.

canopy <- rast("ETH_GlobalCanopyHeight_10m_2020_S27W057_Map.tif")

deep_zoom_ext <- ext(-54.5, -54.4, -25.75, -25.65)

greens <- colorRampPalette(c("#edf8e9", "#a1d99b", "#31a354", "#001709"))

plot(canopy, ext = deep_zoom_ext, col = greens(100),
     main = "Canopy Height, Iguazu National Park",
     axes = TRUE)

To explore the vertical structure of habitats associated with sleeping sites, we cropped the canopy height raster to a buffered extent that covers all recorded locations (±0.02° around minimum and maximum coordinates). This focused raster allows for better spatial correspondence between sleeping sites and vegetation structure.

We then overlaid all sleeping site locations (n = {nrow(moon)}), color-coded by group identity, on the canopy height map. The green gradient reflects vegetation height, with darker shades representing taller canopy. Each monkey group was assigned a distinctive color (e.g., gun in orange, gue in blue), allowing clear visual differentiation of their spatial patterns.

This visualization helps assess whether different groups use sleeping sites in areas of varying canopy height and can inform hypotheses regarding group-specific strategies for predator avoidance, thermoregulation, or social cohesion.

ext_moon <- ext(min(moon$lon)-0.02, max(moon$lon)+0.02,
                min(moon$lat)-0.02, max(moon$lat)+0.02)

canopy_zoomed <- crop(canopy, ext_moon)


group_cols <- c(
    "gun" = "darkorange3",
    "gue" = "midnightblue",
    "mac" = "goldenrod2",
    "spo" = "darkseagreen4")


moon$group <- as.character(moon$group)


plot(canopy_zoomed, col = greens(100),
     main = "Canopy Height Around Sleeping Sites",
     axes = TRUE)


with(moon, points(lon, lat,
                  pch = 21,
                  bg = group_cols[group], 
                  col = "black",          
                  cex = 1.2))

To quantify the vegetation structure at each sleeping site, we extracted canopy height values directly from the Lang et al. (2022) raster using the GPS coordinates of the sites. The sleeping site data were converted into spatial vectors (SpatVector) and projected to match the coordinate reference system (CRS) of the canopy raster.

This allowed us to assign a canopy height value to every sleeping location based on its exact position. The resulting distribution shows that most sites are associated with mid- to high-canopy areas, suggesting a potential preference for structurally complex or elevated sleeping locations.

The density plot below summarizes the canopy height values at all recorded sites. This forms the basis for comparing real sleeping site locations with randomly generated control points in later analyses.

moon_vect <- vect(moon, geom = c("lon", "lat"), crs = crs(canopy))
moon$canopy_height <- extract(canopy, moon_vect)[,2]  

ggplot(moon, aes(x = canopy_height)) +
  geom_density(fill = "darkseagreen", alpha = 0.6) +
  labs(
    title = "Density of Canopy Height at Sleeping Sites",
    x = "Canopy Height (meters)",
    y = "Density"
  ) +
  theme_minimal(base_size = 14)

Distribution of canopy height across observed sleeping sites. The density plot shows the overall tendency of monkeys to select sleeping locations at specific vertical strata within the forest structure.

3.1 Canopy Height vs. Moonlight: Density Visualization

This figure explores the potential relationship between lunar illumination and the vertical position of monkey sleeping sites. Each point represents a site, with jitter added to improve visibility, and background shading representing a 2D kernel density estimate.

The x-axis shows the fraction of the moon illuminated (from new moon at 0.0 to full moon at 1.0), while the y-axis indicates the canopy height (in meters) at the site. The color gradient (darker = higher density) reveals the most frequent combinations of moonlight and canopy structure.

This visualization helps assess whether sleeping sites under darker nights (low moon fraction) tend to occur in taller vegetation, which could suggest behavioral adaptation for predator avoidance or improved concealment during high-risk conditions.

ggplot(moon, aes(x = moon_fraction, y = canopy_height)) +
  geom_point(
    aes(fill = moon_fraction),
    shape = 21,
    color = "black",
    size = 3,
    stroke = 0.3,
    alpha = 0.85) +
  scale_fill_viridis_c(
    option = "B",
    name = "Moon Illum.",
    guide = guide_colorbar(barheight = 8),
    limits = c(0, 1),
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    labels = c("Dark (0.00)", "Low", "Half (0.50)", "High", "Full (1.00)")
  ) +
  coord_cartesian(ylim = c(0, 40)) +
  labs(
    title = "Canopy Height vs. Lunar Illumination",
    x = "Moon Illumination Fraction",
    y = "Canopy Height (m)"
  ) +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.title = element_text(),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 11),
    legend.position = "right"
  )

4 Generating Random Control Points Around Sleeping Sites

To assess whether monkeys are actively selecting taller canopy areas for sleeping under varying moonlight conditions, we generated random control points around each observed sleeping location. For every real site, we simulated 10 random points positioned within a radius of 10 to 100 meters, ensuring they remained within the monkeys’ likely perceptual or movement range for choosing a sleeping location.

We used the {geosphere} package to calculate these new coordinates based on random bearings (0–360°) and distances. The destPoint() function allows us to project a new coordinate a given distance and direction from a known point, taking into account the Earth’s curvature, a key advantage over simple Euclidean methods, especially in geospatial analysis.

These pseudo-absence points share the same date and lunar attributes as the original site, allowing for controlled comparisons while isolating the role of vegetation structure in sleeping site selection.

n_fake <- 10

set.seed(42)
pseudo_data <- moon %>%
  rowwise() %>%
  do({
   
    lat <- .$lat
    lon <- .$lon

    
    distances <- runif(n_fake, min = 40, max = 100)

    
    bearings <- runif(n_fake, min = 0, max = 360)

    
    new_coords <- destPoint(p = c(lon, lat),    
                            b = bearings,       
                            d = distances)      

    tibble(
      lon = new_coords[,1],
      lat = new_coords[,2],
      fecha = .$fecha,
      group = .$group,
      moon_fraction = .$moon_fraction,
      moon_phase_angle = .$moon_phase_angle,
      moon_phase_name = .$moon_phase_name,
      moon_visible = .$moon_visible,
      used = FALSE
    )
  }) %>%
  ungroup()

After generating pseudo-absence (control) points around each real sleeping site, we merged both datasets into a single unified table. The original sleeping sites were labeled with used = TRUE, while the randomly generated points were marked as used = FALSE.

Each entry, whether real or simulated, includes: - Spatial coordinates (lon, lat) - Group identity - Date of observation - Moon illumination metrics (moon_fraction, moon_phase_angle, etc.) - Whether the moon was visible that night - A usage indicator (used)

This combined dataset (moon_all) allows for controlled comparisons between actual monkey sleeping locations and randomly available but unused alternatives. These comparisons are central to testing whether canopy structure or lunar conditions significantly influence sleeping site choice.

moon_used <- moon %>%
  mutate(used = TRUE) %>%
  select(lon, lat, fecha, group, moon_fraction, moon_phase_angle,
         moon_phase_name, moon_visible, used)


moon_all <- bind_rows(moon_used, pseudo_data)

To enable grouped comparisons between used and available points for each sleeping event, we created a unique identifier (id) by combining monkey group and date (group_fecha). This identifier represents a single nightly sleeping site selection event for a given group.

Using dplyr::cur_group_id(), we assigned a numeric ID to each group-date combination, which is shared across the real sleeping site and its 10 associated random control points. This structure is essential for conditional logistic regression or matched comparisons, ensuring that each set of observations is analyzed within its ecological context.

moon_all <- moon_all %>%
  mutate(event_key = paste(group, fecha, sep = "_")) %>%
  group_by(event_key) %>%
  mutate(id = cur_group_id()) %>%
  ungroup() %>%
  select(-event_key)

4.1 Example of a Sleeping Site Selection Event

To illustrate how random control points are spatially distributed around an actual sleeping site, we visualized a single selection event (one group on one night). The focal site (orange) and its 10 random alternatives (gray) are plotted on top of the canopy height raster, cropped to a ~260 m × 260 m window centered on the true location.

A dashed red circle indicates the 100-meter radius within which the control points were generated, while a horizontal arrow provides a visual scale reference (“100 m”). A north indicator is also included for orientation.

This visualization helps communicate the spatial logic behind the pseudo-absence design: real choices are compared to nearby alternatives that share the same lunar and temporal context, isolating the potential effect of vegetation structure.

target_id <- moon_all %>% filter(used) %>% slice(1) %>% pull(id)
subset_points <- moon_all %>% filter(id == target_id)

pt_real <- subset_points %>% filter(used) %>%
  select(lon, lat) %>% unlist()

zoom_ext <- ext(pt_real["lon"] - 0.0013, pt_real["lon"] + 0.0013,
                pt_real["lat"] - 0.0013, pt_real["lat"] + 0.0013)

canopy_focus <- crop(canopy, zoom_ext)

plot(canopy_focus, col = greens(100),
     main = paste("Zoomed Canopy Height - ID", target_id),
     axes = TRUE)

with(subset_points %>% filter(!used),
     points(lon, lat, pch = 21, bg = "gray80", col = "black", cex = 1.4))

with(subset_points %>% filter(used),
     points(lon, lat, pch = 21, bg = "darkorange", col = "black", cex = 1.8))

circle_100 <- destPoint(pt_real, b = seq(0, 360, length.out = 100), d = 100)
lines(circle_100, col = "darkred", lty = 3, lwd = 2)

scale_line <- destPoint(pt_real, b = 90, d = 100) 


x0 <- pt_real["lon"] 
y0 <- pt_real["lat"] + 0.0011 
x1 <- scale_line[1]
y1 <- y0

segments(x0, y0, x1, y1, col = "black", lwd = 2)

arrows(x0, y0, x1, y1, length = 0.05, angle = 90, code = 3, col = "black")

text(x = mean(c(x0, x1)), y = y0 - 0.0001, labels = "100 m", cex = 0.9)

text(x = zoom_ext[1] + 0.0002, y = zoom_ext[4] - 0.0002, labels = "↑ N", cex = 1.1)

Figure: Example of a sleeping site selection event. The real sleeping site (orange) is surrounded by 10 randomly generated control points (gray), constrained within a 100-meter radius (red dashed circle). Background raster shows canopy height. A horizontal arrow denotes a 100 m scale bar, and a north indicator provides spatial orientation.

4.2 Contrasting Sleeping Site Events: Dark vs. Bright Moonlight

To visually compare sleeping site selection patterns under contrasting lunar conditions, we randomly selected 8 sleeping events: 4 occurred under low moon illumination (< 0.2 moon fraction) and 4 under high moon illumination (> 0.8). For each event, the real sleeping site (orange) and its associated pseudo-absence points (gray) are plotted within a 100-meter radius, overlaid on high-resolution canopy height data.

This side-by-side layout enables visual inspection of whether monkeys tend to select taller canopy during darker nights, possibly for predator avoidance or concealment. The consistent structure across panels, including a 100 m scale bar and north orientation, facilitates direct ecological comparison between conditions.

These visualizations set the stage for formal statistical testing of habitat selection as a function of moonlight and vegetation structure.

set.seed(222)

ids_dark <- moon_all %>%
  filter(used, moon_fraction < 0.2) %>%
  distinct(id) %>%
  slice_sample(n = 4) %>%
  pull(id)

ids_bright <- moon_all %>%
  filter(used, moon_fraction > 0.8) %>%
  distinct(id) %>%
  slice_sample(n = 4) %>%
  pull(id)

ids_selected <- c(ids_dark, ids_bright)
breaks_fixed <- seq(10, 26, by = 2)
cols_fixed <- greens(length(breaks_fixed) - 1)
par(mfrow = c(2, 4), mar = c(3, 3, 2, 1), cex = 1.2)

for (target_id in ids_selected) {
  
  subset_points <- moon_all %>% filter(id == target_id)
  pt_real <- subset_points %>% filter(used) %>%
    select(lon, lat) %>% unlist()
  
  zoom_ext <- ext(pt_real["lon"] - 0.0013, pt_real["lon"] + 0.0013,
                  pt_real["lat"] - 0.0013, pt_real["lat"] + 0.0013)
  
  canopy_focus <- crop(canopy, zoom_ext)
  
  plot(canopy_focus,
     col = cols_fixed,
     breaks = breaks_fixed,
     main = paste0("ID ", target_id,
                   "\nMoon: ", round(unique(subset_points$moon_fraction), 2)),
     axes = FALSE,
     legend = TRUE)
  
  with(subset_points %>% filter(!used),
       points(lon, lat, pch = 21, bg = "gray80", col = "black", cex = 1.6))
  
  with(subset_points %>% filter(used),
       points(lon, lat, pch = 21, bg = "darkorange", col = "black", cex = 1.6))
  
  circle_100 <- destPoint(pt_real, b = seq(0, 360, length.out = 100), d = 100)
  lines(circle_100, col = "darkred", lty = 3, lwd = 4)
  
  x0 <- pt_real["lon"] + 0.0003
  y0 <- pt_real["lat"] + 0.0012
  scale_line <- destPoint(c(x0, y0), b = 90, d = 100)
  x1 <- scale_line[1]
  y1 <- y0
  
  segments(x0, y0, x1, y1, col = "black", lwd = 2)
  arrows(x0, y0, x1, y1, length = 0.05, angle = 90, code = 3, col = "black")
  text(mean(c(x0, x1)), y0 - 0.0001, "100 m")
  
  text(x = zoom_ext[1] + 0.0002, y = zoom_ext[4] - 0.0002, labels = "↑ N")
}

Examples of sleeping site selection events under contrasting moonlight conditions. Each panel shows a real sleeping site (orange) and 10 random control points (gray) within a 100 m radius (red dashed circle), overlaid on local canopy height. Four events occurred during low moon illumination (< 0.2), and four during high illumination (> 0.8). A 100-meter scale bar and north indicator are included for spatial reference.

5 Effective Moonlight Calculation

To better capture the ecological relevance of nocturnal lunar illumination at sleeping sites, we computed an index of effective moonlight, combining two key astronomical components:

• Moon illumination fraction: the proportion of the moon’s disc that is lit.
• Moon altitude: its elevation above the horizon, measured at a standardized hour… in this case we used three moments during the night.

We used the suncalc package to calculate moon altitude for each observation, and then estimated effective moonlight as:

\[ \text{moon light effective} = \text{moon_fraction} \times max(0, \sin(\text{moon_altitude})) \]

This formulation ensures that only moons above the horizon contribute to illumination (values below 0 are truncated), and that altitude amplifies the effect of the moon’s light, since a higher moon faces less atmospheric scattering and attenuation.

From an ecological perspective, this metric serves as a biologically meaningful proxy for perceived lunar brightness, which may influence behavioral decisions like sleeping site selection, foraging activity, or predator avoidance among nocturnal or crepuscular species.

To maintain consistency across observations and ensure comparability, we calculate the effective moonlight at three fixed time 20:00 PM (of the previous day) UTC, 00:00 AM and 04:00 AM—each night. This standardized time represents a midpoint when nocturnal primates are typically active and when the moon, if present, is often near its zenith. Using a single reference hour simplifies comparisons by avoiding variation due to changing night lengths, moon rise/set times, or unequal sampling times across dates.

Although this approach does not capture the complete nightly dynamics of lunar brightness, it effectively reflects relative nighttime light conditions across the lunar cycle and seasons. The moon’s altitude—and thus light intensity—varies seasonally (with higher maximum elevations in winter) and cyclically (across lunar phases), so comparing a consistent hour still captures biologically relevant variation at international and annual scales.

Several ecological studies have employed similar approximations. For example, research on black swifts and small songbirds found that both the illuminated fraction and altitude of the moon are strong predictors of nocturnal behavior (e.g., Prinz et al. 2025. In predator–prey systems, models using moon phase and elevation have illuminated habitat selection shifts in felids like ocelots and bobcats (Sergeyev et al. 2023).

We calculated moon altitude at a fixed time (~ 02:00 AM UTC) to standardize our measure of nocturnal lunar illumination. This approach helps us compare different nights and sites under consistent conditions—critical when moonrise, moonset, and night length vary seasonally and by phase.

While this does not capture the full nightly variation in moonlight, it effectively reflects relative illumination across the lunar cycle and seasonal changes in the moon’s elevation. For instance, Śmielak (2023) demonstrated that models combining moon phase and altitude explain > 92% of the variability in actual ground-level moonlight, outperforming phase alone (∼ 60%)

That said, some studies use moon altitude and phase dynamically throughout the night. To better approximate overall lunar illumination, one could sample multiple time points (e.g. 20:00, 00:00, and 04:00 local time) or compute a nightly average of moon_altitude and illumination.

This alternative method—sampling at multiple night-time intervals, would provide a more complete picture of illumination dynamics, especially in studies interested in behavior spanning the entire night. However, it also increases computational complexity and may require integrating over longer periods.

moon_all <- moon_all %>%
  mutate(
    fecha_a = as.POSIXct(paste(fecha - 1, "20:00:00"), tz = "UTC"),
    fecha_b = as.POSIXct(paste(fecha, "00:00:00"), tz = "UTC"),
    fecha_c = as.POSIXct(paste(fecha, "04:00:00"), tz = "UTC"))

pos_a <- getMoonPosition(data = moon_all %>% rename(date = fecha_a) %>% select(date, lat, lon))
pos_b <- getMoonPosition(data = moon_all %>% rename(date = fecha_b) %>% select(date, lat, lon))
pos_c <- getMoonPosition(data = moon_all %>% rename(date = fecha_c) %>% select(date, lat, lon))

alt_scaled_a <- pmax(0, sin(pos_a$altitude))
alt_scaled_b <- pmax(0, sin(pos_b$altitude))
alt_scaled_c <- pmax(0, sin(pos_c$altitude))

moon_all <- moon_all %>%
  mutate(
    moon_altitude_avg = (alt_scaled_a + alt_scaled_b + alt_scaled_c) / 3,
    moon_light_effective_avg = moon_fraction * moon_altitude_avg)

First lets ckeck for a simpler moonlight metric: moon fraction, then the effective moonlight

ggplot(moon_all, aes(x = moon_fraction, y = canopy_height)) +
  geom_jitter(
    aes(fill = used),
    shape = 21, size = 3, stroke = 0.3, 
    alpha = 0.05, color = "black",
    width = 0.03,
    show.legend = TRUE) +
  geom_smooth(
    aes(color = used),
    method = "lm",
    se = TRUE,
    size = 1.9
  ) +
  scale_fill_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c( "Available", "Used"),
    name = "Points"
  ) +
  scale_color_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c( "Available", "Used"),
    name = "Trend"
  ) +
  guides(
    fill = guide_legend(override.aes = list(alpha = 1)),
    color = guide_legend(override.aes = list(size = 1.5))
  ) +
  labs(
    x = "Moon Fraction",
    y = "Canopy Height (m)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(),
    legend.title = element_text(face = "bold"),
    legend.position = "right"
  )

ggplot(moon_all, aes(x = moon_light_effective_avg, y = canopy_height)) +
  geom_jitter(
    aes(fill = used),
    shape = 21, size = 3, stroke = 0.3, 
    alpha = 0.05, color = "black",
    width = 0.03,
    show.legend = TRUE) +
  geom_smooth(
    aes(color = used),
    method = "lm",
    se = TRUE,
    size = 1.9
  ) +
  scale_fill_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c( "Available", "Used"),
    name = "Points"
  ) +
  scale_color_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c( "Available", "Used"),
    name = "Trend"
  ) +
  guides(
    fill = guide_legend(override.aes = list(alpha = 1)),
    color = guide_legend(override.aes = list(size = 1.5))
  ) +
  labs(
    title = "Effective Moonlight vs. Canopy Height",
    x = "Effective Moonlight Index",
    y = "Canopy Height (m)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(),
    legend.title = element_text(face = "bold"),
    legend.position = "right"
  )

Relationship between effective moonlight and canopy height across sleeping site observations. Points represent both used (orange) and available (blue-gray) locations, while lines show linear trends per group. The effective moonlight index combines lunar illumination fraction and altitude to approximate actual night-time light availability.

ggplot(moon_all %>% filter(used == TRUE), aes(x = moon_light_effective_avg, y = canopy_height)) +
  geom_point(
    aes(fill = used),
    shape = 21, size = 5, stroke = 0.3,
    alpha = 0.1, color = "goldenrod4",
    show.legend = TRUE
  ) +
  geom_smooth(
    aes(color = used),
    method = "lm",
    se = TRUE,
    size = 2.9
  ) +
  scale_fill_manual(
    values = c("TRUE" = "#E1B000"),
    labels = c("Used"),
    name = "Points"
  ) +
  scale_color_manual(
    values = c("TRUE" = "goldenrod4"),
    labels = c("Used"),
    name = "Trend"
  ) +
  guides(
    fill = guide_legend(override.aes = list(alpha = 1)),
    color = guide_legend(override.aes = list(size = 1.5))
  ) +
  coord_cartesian(ylim = c(0, 30)) +
  labs(
    title = "Effective Moonlight vs. Canopy Height",
    x = "Effective Moonlight Index",
    y = "Canopy Height (m)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(),
    legend.title = element_text(face = "bold"),
    legend.position = "right"
  )

6 Precipitacion

To incorporate environmental context into our analysis of sleeping site selection, we obtained daily weather data from a meteorological station operated by the Argentine National Meteorological Service (SMN). This station is located approximately 7 km from the study area and provides reliable long-term climatic records.

The dataset includes maximum, minimum, and mean daily temperatures, as well as daily precipitation (mm). We parsed the date field using the lubridate::dmy() function and converted all numeric fields to proper numeric types for further analysis. These data allow us to explore whether rainfall patterns influence the characteristics of sleeping site selection, particularly with regard to vegetation structure and canopy cover.

As precipitation may affect microhabitat conditions, such as humidity, substrate moisture, and thermal regulation, it is plausible that primates modulate their choice of sleeping sites in response to recent rainfall events.

weather <- read.csv("clima.csv")

weather <- weather %>%
  mutate(
    fecha = dmy(fecha),  # formato "d/m/Y"
    t_max = as.numeric(t_max),
    t_min = as.numeric(t_min),
    t_mean = as.numeric(t_mean),
    pp = as.numeric(pp)
  )

6.1 Precipitation Effects on Sleeping Site Structure

To investigate how recent rainfall may influence sleeping site structure, we joined the weather dataset to our spatial records using the fecha (date) field as a key. Specifically, we added the daily precipitation amount (pp, in mm) to the moon_all dataset:

moon_all <- moon_all %>%
  left_join(weather %>% select(fecha, pp), by = "fecha")

We then categorized daily rainfall into three meaningful bins: no rain (0 mm), light rain (1–5 mm), and heavy rain (>5 mm). This allowed us to visually explore whether rainfall levels interact with canopy height and effective moonlight in driving sleeping site selection.

The resulting plot combines:

• Jittered points, filled by rainfall category, to represent variation in conditions across nights
• Separate linear trends for used and available sites, helping assess how rainfall may modulate site choice in different light environments

Together, this visualization provides insight into whether monkeys prefer different canopy structures under varying light and moisture conditions, potentially reflecting thermoregulatory or safety tradeoffs.

moon_all <- moon_all %>%
  left_join(weather %>% select(fecha, pp), by = "fecha")

moon_all <- moon_all %>%
  mutate(pp_cat = cut(
    pp,
    breaks = c(-0.01, 0, 5, Inf),
    labels = c("0 mm", "1–5 mm", ">5 mm")
  ))

ggplot(moon_all, aes(x = moon_light_effective_avg, y = canopy_height)) +
  geom_jitter(
    aes(fill = pp_cat),
    shape = 21, size = 3, stroke = 0.3,
    alpha = 0.5, color = "black", width = 0.03
  ) +
  geom_smooth(
    aes(color = used),
    method = "lm", se = TRUE, size = 1.6
  ) +
  scale_fill_manual(
    values = c("0 mm" = "gray80", "1–5 mm" = "#6baed6", ">5 mm" = "steelblue4"),
    name = "Precipitation"
  ) +
  scale_color_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c("Used", "Available"),
    name = "Trend"
  ) +
  labs(
    title = "Canopy Height vs. Effective Moonlight and Precipitation",
    x = "Effective Moonlight Index",
    y = "Canopy Height (m)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    legend.title = element_text(face = "bold")
  )

Figure: Relationship between canopy height and effective moonlight under varying precipitation levels. Points are colored by precipitation categories recorded on the same date: no rain (gray), light rain (blue), and moderate to heavy rain (dark blue). Trend lines show linear patterns for used (yellow) and available (blue) locations. While moonlight may influence site selection, canopy height patterns do not appear strongly driven by precipitation.

ggplot(moon_all, aes(x = pp, y = canopy_height, color = used)) +
  geom_jitter(
    shape = 21, fill = "#6baed6", color = "black",
    alpha = 0.1, width = 0.3, size = 2.6
  ) +
  geom_smooth(method = "lm", se = TRUE, size = 2.5) +
  scale_color_manual(
    values = c("TRUE" = "#E1B000", "FALSE" = "#1F3552"),
    labels = c("Used", "Available"),
    name = "Trend"
  ) +
  labs(
    title = "Canopy Height vs. Precipitation",
    x = "Precipitation (mm)",
    y = "Canopy Height (m)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.title = element_text(face = "bold")
  )

Figure: Relationship between daily precipitation and canopy height at sleeping sites. Each point represents a used (yellow) or available (blue) location on a given date, with precipitation recorded at a nearby weather station. Smoother lines represent linear trends for each category. Although precipitation varies substantially, no strong association with canopy height is evident across categories.