Load packages

library(here)
library(readr)
library(dplyr)
library(ggplot2)
library(readxl)

Load data

Path is built from the directory of this Rmd (when run in RStudio) so it works even if another project is open. Otherwise falls back to here() (project root). Same layout works for everyone and after cloning from GitHub.

# Use this script's directory so the correct data is found regardless of open project
script_dir <- if (requireNamespace("rstudioapi", quietly = TRUE)) {
  tryCatch(
    dirname(rstudioapi::getActiveDocumentContext()$path),
    error = function(e) here()
  )
} else {
  here()
}
path_data <- file.path(script_dir, "data_processed_environment_nms_20250603.csv")
if (!file.exists(path_data)) stop("Data file not found at: ", path_data, "\nEnsure this Rmd and the CSV are in the same folder (e.g. Grassworks_Data).")
sites <- read_csv(path_data, show_col_types = FALSE)
head(sites)

Filter sites

sites <- sites %>%
  filter(region == "south", !is.na(rest.age))
head(sites)
length(unique(sites$id.site)) 
n_per_site <- sites %>% count(id.site, name = "n_obs")

Part 2: Species data

Load abundance data and keep only rows for sites that passed the site filter (region south, rest.age not NA).

path_species <- file.path(script_dir, "data_processed_abundances_t_20250306.csv")
if (!file.exists(path_species)) stop("Species file not found at: ", path_species)
species <- read_csv(path_species, show_col_types = FALSE)
# Restrict to the same id.site as in filtered sites
sites_kept <- unique(sites$id.site)
species <- species %>% filter(id.site %in% sites_kept)
head(species)

Species list: seed_to_cover

Load the Excel sheet and build a df of species sown in greenhous (with seed_to_cover == 1).

# Path: one folder up from script dir (Merlin_Ruehrer), same for all machines
path_seed_to_cover <- file.path(script_dir, "..", "Arten_Einzelartenexperiment_seed_to_cover.xlsx")
if (!file.exists(path_seed_to_cover)) stop("Excel file not found at: ", path_seed_to_cover)
seed_to_cover_raw <- read_excel(path_seed_to_cover)
specs_cover <- subset(seed_to_cover_raw, seed_to_cover == 1)
nrow(specs_cover)
## [1] 82

Species frequency in species (by plots and by sites)

For each species (name.plant) in df species:

n_plots_total <- n_distinct(species$id.plot)
n_sites_total <- n_distinct(species$id.site)

species_freq_all <- species %>%
  group_by(name.plant) %>%
  summarise(
    n_plots = n(),                                    # all rows = plot-level occurrences
    n_sites = n_distinct(id.site),                    # sites where species occurs (≥1 row per site)
    .groups = "drop"
  ) %>%
  mutate(
    freq_plots = n_plots / n_plots_total,
    freq_sites = n_sites / n_sites_total
  )

View(species_freq_all)

Add the frequencies and name.plant to specs_cover, matching by Artname (= name.plant in species data).
An exact join misses cases where one side has infraspecific notation (e.g. Dactylis glomerata s. str.) and the other does not. We use a normalized name (genus + species, stripping “s. str.”, “s. l.”, “agg.”, etc.) so that both match.

# Normalize plant name to genus + species (strip infraspecific / aggregate notation)
normalize_plant_name <- function(x) {
  x %>%
    stringr::str_remove_all("\\s+s\\.\\s*str\\.?|\\s+s\\.\\s*l\\.?|\\s+agg\\.?|\\s+subsp\\.?|\\s+var\\.?") %>%
    stringr::str_squish() %>%
    # Keep first two words (genus species) as canonical match key
    stringr::str_extract("^\\S+\\s+\\S+") %>%
    trimws()
}

# Build species frequencies by normalized name (sum if multiple name.plant map to same binomial)
species_freq_norm <- species_freq_all %>%
  mutate(name_norm = normalize_plant_name(name.plant)) %>%
  filter(!is.na(name_norm)) %>%
  group_by(name_norm) %>%
  summarise(
    n_plots_norm = sum(n_plots),
    n_sites_norm = sum(n_sites),
    freq_plots_norm = n_plots_norm / n_plots_total,
    freq_sites_norm = n_sites_norm / n_sites_total,
    .groups = "drop"
  )

# Exact-match columns under fixed names so coalesce always has both sides
species_freq_exact <- species_freq_all %>%
  rename(
    n_plots_exact = n_plots,
    n_sites_exact = n_sites,
    freq_plots_exact = freq_plots,
    freq_sites_exact = freq_sites
  )

# Join: normalized first (always get _norm), then exact (get _exact); coalesce to final cols
specs_cover <- specs_cover %>%
  mutate(Artname_norm = normalize_plant_name(Artname)) %>%
  left_join(species_freq_norm, by = c("Artname_norm" = "name_norm")) %>%
  left_join(
    species_freq_exact,
    by = c("Artname" = "name.plant")
  ) %>%
  mutate(
    n_plots = coalesce(n_plots_exact, n_plots_norm),
    n_sites = coalesce(n_sites_exact, n_sites_norm),
    freq_plots = coalesce(freq_plots_exact, freq_plots_norm),
    freq_sites = coalesce(freq_sites_exact, freq_sites_norm),
    name.plant = Artname
  ) %>%
  select(-Artname_norm, -n_plots_norm, -n_sites_norm, -freq_plots_norm, -freq_sites_norm,
         -n_plots_exact, -n_sites_exact, -freq_plots_exact, -freq_sites_exact)

View(specs_cover)

I want to keep those species that occur in at least 10 sites

specs_cover_10 <- subset(specs_cover, n_sites >= 10)
length(unique(specs_cover_10$name.plant))
## [1] 44
unique(specs_cover_10$name.plant)
##  [1] "Achillea millefolium"    "Agrimonia eupatoria"    
##  [3] "Agrostis capillaris"     "Ajuga reptans"          
##  [5] "Alopecurus pratensis"    "Anthoxanthum odoratum"  
##  [7] "Arrhenatherum elatius"   "Briza media"            
##  [9] "Campanula patula"        "Carex flacca"           
## [11] "Carum carvi"             "Crepis biennis"         
## [13] "Cynosurus cristatus"     "Dactylis glomerata"     
## [15] "Daucus carota"           "Festuca pratensis"      
## [17] "Festuca rubra"           "Galium verum"           
## [19] "Holcus lanatus"          "Hypericum perforatum"   
## [21] "Hypochaeris radicata"    "Knautia arvensis"       
## [23] "Lathyrus pratensis"      "Leontodon hispidus"     
## [25] "Lotus corniculatus"      "Medicago lupulina"      
## [27] "Phleum pratense"         "Plantago lanceolata"    
## [29] "Poa pratensis"           "Prunella vulgaris"      
## [31] "Ranunculus acris"        "Ranunculus bulbosus"    
## [33] "Rhinanthus minor"        "Rumex acetosa"          
## [35] "Salvia pratensis"        "Sanguisorba minor"      
## [37] "Sanguisorba officinalis" "Silene vulgaris"        
## [39] "Stellaria graminea"      "Tragopogon pratensis"   
## [41] "Trifolium pratense"      "Trisetum flavescens"    
## [43] "Veronica chamaedrys"     "Vicia cracca"

44 Arten kommen auf mindestens 10 Flächen aus dem Grassworks Projekt vor (nur südliche Flächen mit Informationen zum Alter)

family,n - Poaceae,13 - Asteraceae,5 - Fabaceae,5 - Lamiaceae,3 - Rosaceae,3 - Apiaceae,2 - Caryophyllaceae,2 - Plantaginaceae,2 - Ranunculaceae,2 - Campanulaceae,1 - Caprifoliaceae,1 - Cyperaceae,1 - Hypericaceae,1 - Orobanchaceae,1 - Polygonaceae,1 - Rubiaceae,1

View(specs_cover_10)

in addition: those that are seeded at Inn Dikes to have the opportunity to extend the ages of the grasslands to 1 year (gw-data: 2-36 yrs)

specs_cover_11 <- specs_cover_10 %>%
  filter(InnDeiche == 1)
unique(specs_cover_11$name.plant)
##  [1] "Achillea millefolium"    "Agrimonia eupatoria"    
##  [3] "Anthoxanthum odoratum"   "Arrhenatherum elatius"  
##  [5] "Briza media"             "Campanula patula"       
##  [7] "Carum carvi"             "Crepis biennis"         
##  [9] "Dactylis glomerata"      "Festuca rubra"          
## [11] "Galium verum"            "Holcus lanatus"         
## [13] "Hypericum perforatum"    "Hypochaeris radicata"   
## [15] "Knautia arvensis"        "Lathyrus pratensis"     
## [17] "Leontodon hispidus"      "Lotus corniculatus"     
## [19] "Medicago lupulina"       "Phleum pratense"        
## [21] "Plantago lanceolata"     "Poa pratensis"          
## [23] "Prunella vulgaris"       "Ranunculus acris"       
## [25] "Rhinanthus minor"        "Rumex acetosa"          
## [27] "Salvia pratensis"        "Sanguisorba minor"      
## [29] "Sanguisorba officinalis" "Silene vulgaris"        
## [31] "Stellaria graminea"      "Tragopogon pratensis"   
## [33] "Trifolium pratense"      "Trisetum flavescens"
length(unique(specs_cover_11$name.plant))
## [1] 34

Step 1: Define the target species list

Use your list as a character vector.

target_species <- specs_cover_11$Artname
length(target_species)  
## [1] 34

Step 2: Match species data to the target list (handle name variants)

In species, the column is name.plant; it may contain infraspecific names (e.g. “Dactylis glomerata s. str.”). We normalize to genus + species so they match the target list.

# Use same normalizer as earlier in this script
species_target <- species %>%
  mutate(name_norm = normalize_plant_name(name.plant)) %>%
  filter(name_norm %in% normalize_plant_name(target_species))

# How many plot × species records (and distinct sites) have at least one target species?
nrow(species_target)
## [1] 1657
n_distinct(species_target$id.site)
## [1] 32

Step 3: Per site – which target species occur there, and how many?

For each id.site, count distinct target species (a species counts once per site even if in several plots).

sites_with_targets <- species_target %>%
  group_by(id.site) %>%
  summarise(
    n_target_species = n_distinct(name_norm),
    target_species_list = list(sort(unique(name_norm))),
    .groups = "drop"
  ) %>%
  arrange(desc(n_target_species))

# Preview: sites with most target species
head(sites_with_targets, 10)
## # A tibble: 10 × 3
##    id.site n_target_species target_species_list
##    <chr>              <int> <list>             
##  1 S_HTH                 30 <chr [30]>         
##  2 S_BUC                 27 <chr [27]>         
##  3 S_HUB                 27 <chr [27]>         
##  4 S_WNK                 27 <chr [27]>         
##  5 S_AUH                 26 <chr [26]>         
##  6 S_OST                 26 <chr [26]>         
##  7 S_KRT                 25 <chr [25]>         
##  8 S_OSW                 25 <chr [25]>         
##  9 S_SBN                 24 <chr [24]>         
## 10 S_TIE                 24 <chr [24]>

Step 4: Site selection – which sites to visit to be most efficient

10 sites per species with as little effort as possible (for all 31 species until now)

# 2) Normalize names in both target list and species df
target_norm <- normalize_plant_name(target_species)

species_target <- species %>%
  mutate(name_norm = normalize_plant_name(name.plant)) %>%
  filter(name_norm %in% target_norm)
# 3) Double check: does each target species occur in ≥ 10 sites ?
species_site_counts <- species_target %>%
  group_by(name_norm) %>%
  summarise(
    n_sites_total = n_distinct(id.site),
    .groups = "drop"
  )

View(species_site_counts)
# 4) For each site: which target species (normalized) occur there?
site_species <- species_target %>%
  distinct(id.site, name_norm)

site_list <- site_species %>%
  group_by(id.site) %>%
  summarise(
    species_here = list(sort(unique(name_norm))),
    .groups = "drop"
  )
# ---------- choose sites so that each species gets 10 sites ----------

required_per_species <- 10  # change this number if needed

# Current count of how many chosen sites each species has
current_counts <- setNames(
  rep(0L, length(target_norm)),
  sort(unique(target_norm))
)

sites_remaining <- site_list
sites_chosen <- character(0)

repeat {
  # How much does each remaining site help towards the 10-site goal?
 gains <- vapply(
  seq_len(nrow(sites_remaining)),
  function(i) {
    sp <- sites_remaining$species_here[[i]]
    # For species at this site, how much do we increase the count,
    # but not beyond "required_per_species"?
    sum(pmax(
      0L,
      pmin(required_per_species, current_counts[sp] + 1L) - current_counts[sp]
    ))
  },
  numeric(1))

  # If no site adds anything, break
  if (all(gains == 0L)) break

  # Choose site with maximum gain
  best_idx <- which.max(gains)
  best_site <- sites_remaining$id.site[best_idx]
  best_species <- sites_remaining$species_here[[best_idx]]

  sites_chosen <- c(sites_chosen, best_site)

  # Update counts: each species at this site gains 1 visit
  for (sp in best_species) {
    current_counts[sp] <- min(required_per_species, current_counts[sp] + 1L)
  }

  # Remove chosen site from remaining
  sites_remaining <- sites_remaining[-best_idx, ]

  # Stop if all species reached 10 sites
  if (all(current_counts >= required_per_species)) break
}
# Sites to visit (order: best first, second best, …)
sites_chosen
##  [1] "S_HTH" "S_BUC" "S_HUB" "S_WNK" "S_AUH" "S_OST" "S_KRT" "S_OSW" "S_SBN"
## [10] "S_TIE" "S_ACR" "S_GIG" "S_GSH" "S_JAU" "S_MBM" "S_FCH" "S_HLZ" "S_NOZ"
## [19] "S_JHB" "S_GTZ" "S_POP" "S_SHH" "S_ECH" "S_GSR"
length(sites_chosen)
## [1] 24
  • if we measure all species that occur in greenhouse and Inn Dikes and in at least 10 sites of gw-project, we have to visit 24 sites maximum plus Inn Dikes.

We want age as predictor (x-variable) - how well is the age gradient (proxy for community completeness, evenness etc) distributed within those 24 sites?

# Join with the sites df to get rest.age
sites_chosen_info <- sites %>%
  filter(id.site %in% sites_chosen) %>%
  distinct(id.site, rest.age)

View(sites_chosen_info)

# Histogram of restoration ages represented by the chosen sites
hist(
  sites_chosen_info$rest.age,
  breaks = 20,                       # adjust if needed
  main   = "Restoration ages of chosen sites",
  xlab   = "rest.age",
  col    = "grey80",
  border = "white"
)

## can we obtain well distributed ages spans for each species??

# Restrict to chosen sites, add rest.age
species_chosen <- species_target %>%
  filter(id.site %in% sites_chosen) %>%
  left_join(
    sites %>% select(id.site, rest.age),
    by = "id.site"
  )

# For each target species: which restoration ages (and how many) do we cover?
species_age_coverage <- species_chosen %>%
  distinct(name_norm, id.site, rest.age) %>%      # 1 row per species × site
  group_by(name_norm) %>%
  summarise(
    ages_covered = list(sort(unique(rest.age))),  # vector of ages per species
    n_ages       = length(unique(rest.age)),      # how many distinct ages
    n_sites      = n_distinct(id.site),           # how many sites overall
    .groups = "drop"
  )

View(species_age_coverage)
# One row per species × site × age
species_age_long <- species_chosen %>%
  distinct(name_norm, id.site, rest.age)

# Faceted histogram: one panel per species
ggplot(species_age_long, aes(x = rest.age)) +
  geom_histogram(binwidth = 5, boundary = 0, color = "white", fill = "grey60") +
  facet_wrap(~ name_norm, scales = "free_y") +
  labs(
    x = "rest.age",
    y = "Number of sites",
    title = "Restoration ages per species in chosen sites"
  ) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Arten mit relativ guter “Altersverteilung” in den Flächen (visuell)

Hauptbluete FloraWeb Monate

earlier:

later:

weitere Überlegungen zu diesen Arten

todo: check wie viele Sites angefahren werden müssten mit der finalen Auswahl an Arten

maximal (hier mit allen Arten): 24 sites (s.o.)

todo: mit der finalen Auswahl an Arten prüfen, auf wie vielen PLOTS (hier erstmal nur sites als Kriterium genutzt) pro site diese auch vorhanden sind.

todo: Testlauf Freiland möglichst bald: Manschette um “Zielart”; Foto; Kiste; anhand dessen entscheiden: Wie viele Arten machbar und wie viele obs. per spec within each plot / site.