From 2bd594c3b9f7a5cd4d28869a151825a58a8771a5 Mon Sep 17 00:00:00 2001
From: chross22 <52218551+chross22@users.noreply.github.com>
Date: Thu, 6 Aug 2026 15:11:39 -0400
Subject: [PATCH] Add the Labrador Current retroflection index, with its
 citation

Jutras et al. publish their index as source data with the paper, so this
fetches the authors' own values rather than reimplementing them. It
joins the existing catalog as LCR alongside NAO, AO, AMO, and PDO.

It is the odd one out among them and often the more directly useful: the
other four are atmospheric patterns, while this describes a current -
how much of the Labrador Current turns east at the Grand Banks instead
of continuing along the shelf, and so how much cold, fresh,
oxygen-rich water reaches the Scotian Shelf and Gulf of Maine. For shelf
water properties that is a shorter causal chain than the NAO.

Being the output of one study rather than an operational product, it
carries a citation that travels with it: a new `reference` field on the
catalog entry, surfaced by index_dictionary() both as a column and in
its printed output, and repeated in the fetch_climate_index() docs and
the README. The NOAA indices have no single paper to point at, so their
reference is NA rather than invented.

Two format details, both established from the file rather than assumed:

  - The dates are decimal years on a fixed 365-day year. The step
    between rows is 1/365 to eleven significant figures; reading them as
    365.25 would walk the derived dates off by several days across the
    record and move values into neighbouring months.
  - Daily values are averaged to monthly, since every other index here
    is monthly and attach_climate_index() joins on year and month. The
    underlying index is already 12-month smoothed, so little is lost.

Months backed by fewer than half their days are dropped. The 365-day
convention rolls the final record into a January of its own, which would
otherwise appear as a confident-looking January 2015 value resting on a
single day - the failure mode min_coverage guards against elsewhere.

Verified against the live file: 260 monthly values, May 1993 to December
2014, attaching correctly to observations.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 R/climate_indices.R                   | 131 +++++++++++++++++++++++++-
 README.Rmd                            |  34 ++++++-
 README.md                             |  35 ++++++-
 man/fetch_climate_index.Rd            |  10 ++
 man/parse_decimal_year_csv.Rd         |  41 ++++++++
 tests/testthat/test-climate-indices.R |  64 +++++++++++++
 6 files changed, 311 insertions(+), 4 deletions(-)
 create mode 100644 man/parse_decimal_year_csv.Rd

diff --git a/R/climate_indices.R b/R/climate_indices.R
index f299888..1f9fddb 100644
--- a/R/climate_indices.R
+++ b/R/climate_indices.R
@@ -53,6 +53,27 @@ climate_indices <- function() {
       description = paste("Leading mode of North Pacific SST variability.",
                           "Included for completeness; of limited relevance to",
                           "Atlantic shelf systems.")
+    ),
+    LCR = list(
+      label = "Labrador Current retroflection",
+      source = "Jutras et al. 2023, Nature Communications",
+      url = paste0("https://static-content.springer.com/esm/",
+                   "art%3A10.1038%2Fs41467-023-38321-y/MediaObjects/",
+                   "41467_2023_38321_MOESM5_ESM.csv"),
+      format = "decimal_year_csv",
+      reference = paste("Jutras M, Dufour CO, Mucci A, Talbot LC (2023)",
+                        "Large-scale control of the retroflection of the",
+                        "Labrador Current. Nature Communications 14:2623.",
+                        "doi:10.1038/s41467-023-38321-y"),
+      description = paste("How much of the Labrador Current turns eastward at",
+                          "the Grand Banks instead of continuing southwest along",
+                          "the shelf. Positive values mean stronger retroflection,",
+                          "so less cold, fresh, oxygen-rich Labrador water reaches",
+                          "the Scotian Shelf and Gulf of Maine. Unlike the",
+                          "atmospheric indices here it describes a current rather",
+                          "than a pressure or temperature pattern, which makes it",
+                          "the more direct predictor of shelf water properties.",
+                          "Covers 1993-2014 only.")
     )
   )
 }
@@ -68,7 +89,12 @@ index_dictionary <- function() {
   dictionary <- do.call(rbind, lapply(names(catalog), function(name) {
     entry <- catalog[[name]]
     data.frame(name = name, label = entry$label, source = entry$source,
-               url = entry$url, description = entry$description,
+               url = entry$url,
+               # Indices published with a paper carry its citation. The
+               # operational ones from NOAA have no single paper to point at,
+               # so this is empty for them rather than invented.
+               reference = entry$reference %||% NA_character_,
+               description = entry$description,
                stringsAsFactors = FALSE)
   }))
   class(dictionary) <- c("datamatch_index_dictionary", "data.frame")
@@ -85,7 +111,16 @@ print.datamatch_index_dictionary <- function(x, ...) {
   cat(strrep("-", 62), "\n", sep = "")
   print(flat[c("name", "label", "source")], row.names = FALSE, right = FALSE)
   cat("\nThese have no spatial dimension: one value per month, basin-wide.\n")
-  cat("Sources: as.data.frame(index_dictionary())$url\n")
+
+  cited <- flat[!is.na(flat$reference), ]
+  if (nrow(cited) > 0) {
+    cat("\nCite when used:\n")
+    for (i in seq_len(nrow(cited))) {
+      cat("  ", cited$name[i], ": ", cited$reference[i], "\n", sep = "")
+    }
+  }
+
+  cat("\nSources: as.data.frame(index_dictionary())$url\n")
   invisible(x)
 }
 
@@ -97,6 +132,18 @@ print.datamatch_index_dictionary <- function(x, ...) {
 #' The published files use fixed-width year-by-month tables with provider-specific
 #' missing-value codes, which are parsed here into one row per month.
 #'
+#' @section Citing an index:
+#' `LCR` is the published output of a specific study rather than an operational
+#' product, and should be cited when used:
+#'
+#' > Jutras M, Dufour CO, Mucci A, Talbot LC (2023) Large-scale control of the
+#' > retroflection of the Labrador Current. *Nature Communications* **14**:2623.
+#' > \doi{10.1038/s41467-023-38321-y}
+#'
+#' The series is the source data published with that paper's Figure 3, fetched
+#' from the journal rather than recomputed, so the values are the authors' own.
+#' `as.data.frame(index_dictionary())$reference` carries this at runtime.
+#'
 #' @param index an index name from [climate_indices()]
 #' @param years years to keep; `NULL` keeps the whole record
 #' @param url override the catalog URL, if a provider has moved the file
@@ -145,6 +192,13 @@ fetch_climate_index <- function(index, years = NULL, url = NULL) {
 #' @return a data frame with `YEAR`, `MONTH`, and the named value column
 #' @keywords internal
 parse_index_table <- function(lines, format, name) {
+  # Not every published index is a year-by-month table. The retroflection index
+  # is a daily series in decimal years, so it is parsed separately rather than
+  # bent into the shape of the CPC and PSL files.
+  if (format == "decimal_year_csv") {
+    return(parse_decimal_year_csv(lines, name))
+  }
+
   lines <- trimws(lines)
   lines <- lines[nzchar(lines)]
 
@@ -196,6 +250,79 @@ parse_index_table <- function(lines, format, name) {
   series
 }
 
+#' Parse a daily index published as decimal years
+#'
+#' The format the Jutras et al. retroflection index is published in: a couple of
+#' figure-caption lines, a `Date,<label>` header, then one row per day with the
+#' date as a decimal year.
+#'
+#' @section Two conversions worth stating:
+#' **Decimal years use a fixed 365-day year here.** The step between consecutive
+#' rows is 1/365 to eleven significant figures, so leap days are not represented
+#' and a date is recovered as day `round(fraction * 365)` of its year. Assuming
+#' 365.25 instead would walk the derived dates off by up to several days over the
+#' record.
+#'
+#' **Daily values are averaged to monthly.** Every other index in this package is
+#' monthly, and [attach_climate_index()] joins on year and month, so a daily
+#' series has nothing to join to. The underlying index is already smoothed with a
+#' 12-month rolling mean, so monthly averaging discards very little.
+#'
+#' Months backed by fewer than half their days are dropped rather than reported.
+#' The record's first and last months are partial by construction, and the fixed
+#' 365-day year rolls the final record into a January of its own — which would
+#' otherwise surface as a confident-looking monthly value resting on a single day.
+#'
+#' @param lines the downloaded file's lines
+#' @param name what to call the value column
+#' @return a data frame with `YEAR`, `MONTH`, and the named value column
+#' @keywords internal
+parse_decimal_year_csv <- function(lines, name) {
+  lines <- trimws(lines)
+  # A data row is two comma-separated numbers. The figure-caption and header
+  # lines fail that, which is steadier than skipping a fixed number of lines.
+  is_data <- grepl("^-?[0-9.]+,-?[0-9.eE+-]+$", lines)
+  rows <- lines[is_data]
+
+  if (length(rows) == 0) {
+    stop("No data rows found in the index file; its format may have changed.",
+         call. = FALSE)
+  }
+
+  fields <- do.call(rbind, strsplit(rows, ","))
+  decimal_year <- as.numeric(fields[, 1])
+  value <- as.numeric(fields[, 2])
+  keep <- !is.na(decimal_year) & !is.na(value)
+  decimal_year <- decimal_year[keep]
+  value <- value[keep]
+
+  year <- floor(decimal_year)
+  date <- as.Date(paste0(year, "-01-01")) + round((decimal_year - year) * 365)
+
+  daily <- data.frame(YEAR = as.integer(format(date, "%Y")),
+                      MONTH = as.integer(format(date, "%m")),
+                      value = value)
+
+  monthly <- stats::aggregate(value ~ YEAR + MONTH, data = daily,
+                              FUN = function(z) c(mean(z), length(z)))
+  monthly <- data.frame(YEAR = monthly$YEAR, MONTH = monthly$MONTH,
+                        value = monthly$value[, 1], days = monthly$value[, 2])
+
+  # A month backed by a handful of days is not that month's value. The record's
+  # first and last months are partial by construction, and the fixed 365-day
+  # year rolls the final record into a January of its own - which would
+  # otherwise appear as a confident-looking value resting on one day. Half a
+  # month is the same bar min_coverage sets by default elsewhere here.
+  full <- as.numeric(lubridate::days_in_month(
+    lubridate::ymd(paste(monthly$YEAR, monthly$MONTH, 1, sep = "-"))))
+  monthly <- monthly[monthly$days >= full / 2, ]
+
+  monthly <- monthly[order(monthly$YEAR, monthly$MONTH), c("YEAR", "MONTH", "value")]
+  names(monthly)[names(monthly) == "value"] <- name
+  rownames(monthly) <- NULL
+  monthly
+}
+
 #' Missing-value code declared in a PSL index file
 #'
 #' PSL files state their own missing-value code on a line after the data, rather
diff --git a/README.Rmd b/README.Rmd
index 0cd604a..5e62448 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -436,7 +436,7 @@ so every observation in a month receives the same number:
 ```{r indices, eval = FALSE}
 observations <- attach_climate_index(observations, c("NAO", "AMO"))
 
-index_dictionary()          # NAO, AO, AMO, PDO
+index_dictionary()          # NAO, AO, AMO, PDO, LCR
 ```
 
 That makes them a different kind of covariate from local temperature. They carry
@@ -445,6 +445,38 @@ within a region conditions are better. Useful for interannual questions ("was
 this a warm-regime year?"), useless for spatial ones: a model given only indices
 cannot produce a map.
 
+#### The Labrador Current retroflection index
+
+`LCR` is the odd one out and often the most directly useful of the five. The
+other four are atmospheric patterns; this one describes a *current* — how much
+of the Labrador Current turns eastward at the Grand Banks instead of continuing
+southwest along the shelf. Positive values mean stronger retroflection, so less
+cold, fresh, oxygen-rich Labrador water reaching the Scotian Shelf and Gulf of
+Maine. For shelf water properties that is a shorter causal chain than the NAO.
+
+```{r lcr, eval = FALSE}
+lcr <- fetch_climate_index("LCR")        # monthly, 1993-2014
+observations <- attach_climate_index(observations, "LCR")
+```
+
+It is the published output of one study rather than an operational product, so
+**cite it when you use it**:
+
+> Jutras M, Dufour CO, Mucci A, Talbot LC (2023) Large-scale control of the
+> retroflection of the Labrador Current. *Nature Communications* **14**:2623.
+> <https://doi.org/10.1038/s41467-023-38321-y>
+
+`index_dictionary()` prints that citation, and
+`as.data.frame(index_dictionary())$reference` carries it at runtime.
+
+Two limits worth knowing. It **covers 1993–2014 only**, so it cannot be attached
+to recent observations. And the series is fetched from the paper's published
+source data rather than recomputed — the authors derived it by tracking virtual
+particles through GLORYS12V1 with
+[OceanParcels](https://oceanparcels.org/), which is a Lagrangian modelling
+exercise rather than something this package reproduces. Extending the record
+past 2014 means redoing that computation, not calling a different function.
+
 One interaction to know: `attach_climate_index()` joins on year and month, and
 `upscale_time(to = "year")` stamps its output `MONTH = 1`. Attaching an index to
 annual data would therefore give every year January's value. Aggregate the index
diff --git a/README.md b/README.md
index df22606..4658457 100644
--- a/README.md
+++ b/README.md
@@ -508,7 +508,7 @@ so every observation in a month receives the same number:
 ``` r
 observations <- attach_climate_index(observations, c("NAO", "AMO"))
 
-index_dictionary()          # NAO, AO, AMO, PDO
+index_dictionary()          # NAO, AO, AMO, PDO, LCR
 ```
 
 That makes them a different kind of covariate from local temperature. They carry
@@ -517,6 +517,39 @@ within a region conditions are better. Useful for interannual questions ("was
 this a warm-regime year?"), useless for spatial ones: a model given only indices
 cannot produce a map.
 
+#### The Labrador Current retroflection index
+
+`LCR` is the odd one out and often the most directly useful of the five. The
+other four are atmospheric patterns; this one describes a *current* — how much
+of the Labrador Current turns eastward at the Grand Banks instead of continuing
+southwest along the shelf. Positive values mean stronger retroflection, so less
+cold, fresh, oxygen-rich Labrador water reaching the Scotian Shelf and Gulf of
+Maine. For shelf water properties that is a shorter causal chain than the NAO.
+
+
+``` r
+lcr <- fetch_climate_index("LCR")        # monthly, 1993-2014
+observations <- attach_climate_index(observations, "LCR")
+```
+
+It is the published output of one study rather than an operational product, so
+**cite it when you use it**:
+
+> Jutras M, Dufour CO, Mucci A, Talbot LC (2023) Large-scale control of the
+> retroflection of the Labrador Current. *Nature Communications* **14**:2623.
+> <https://doi.org/10.1038/s41467-023-38321-y>
+
+`index_dictionary()` prints that citation, and
+`as.data.frame(index_dictionary())$reference` carries it at runtime.
+
+Two limits worth knowing. It **covers 1993–2014 only**, so it cannot be attached
+to recent observations. And the series is fetched from the paper's published
+source data rather than recomputed — the authors derived it by tracking virtual
+particles through GLORYS12V1 with
+[OceanParcels](https://oceanparcels.org/), which is a Lagrangian modelling
+exercise rather than something this package reproduces. Extending the record
+past 2014 means redoing that computation, not calling a different function.
+
 One interaction to know: `attach_climate_index()` joins on year and month, and
 `upscale_time(to = "year")` stamps its output `MONTH = 1`. Attaching an index to
 annual data would therefore give every year January's value. Aggregate the index
diff --git a/man/fetch_climate_index.Rd b/man/fetch_climate_index.Rd
index 44e3b5f..c3e7a0a 100644
--- a/man/fetch_climate_index.Rd
+++ b/man/fetch_climate_index.Rd
@@ -24,6 +24,16 @@ series.
 The published files use fixed-width year-by-month tables with provider-specific
 missing-value codes, which are parsed here into one row per month.
 }
+\section{Citing an index}{
+
+\code{LCR} is the published output of a specific study rather than an operational
+product, and should be cited when used:Jutras M, Dufour CO, Mucci A, Talbot LC (2023) Large-scale control of theretroflection of the Labrador Current. Nature Communications 14:2623.\doi{10.1038/s41467-023-38321-y}
+
+The series is the source data published with that paper's Figure 3, fetched
+from the journal rather than recomputed, so the values are the authors' own.
+\code{as.data.frame(index_dictionary())$reference} carries this at runtime.
+}
+
 \examples{
 \dontrun{
 nao <- fetch_climate_index("NAO", years = 2000:2020)
diff --git a/man/parse_decimal_year_csv.Rd b/man/parse_decimal_year_csv.Rd
new file mode 100644
index 0000000..0e8bd96
--- /dev/null
+++ b/man/parse_decimal_year_csv.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/climate_indices.R
+\name{parse_decimal_year_csv}
+\alias{parse_decimal_year_csv}
+\title{Parse a daily index published as decimal years}
+\usage{
+parse_decimal_year_csv(lines, name)
+}
+\arguments{
+\item{lines}{the downloaded file's lines}
+
+\item{name}{what to call the value column}
+}
+\value{
+a data frame with \code{YEAR}, \code{MONTH}, and the named value column
+}
+\description{
+The format the Jutras et al. retroflection index is published in: a couple of
+figure-caption lines, a \verb{Date,<label>} header, then one row per day with the
+date as a decimal year.
+}
+\section{Two conversions worth stating}{
+
+\strong{Decimal years use a fixed 365-day year here.} The step between consecutive
+rows is 1/365 to eleven significant figures, so leap days are not represented
+and a date is recovered as day \code{round(fraction * 365)} of its year. Assuming
+365.25 instead would walk the derived dates off by up to several days over the
+record.
+
+\strong{Daily values are averaged to monthly.} Every other index in this package is
+monthly, and \code{\link[=attach_climate_index]{attach_climate_index()}} joins on year and month, so a daily
+series has nothing to join to. The underlying index is already smoothed with a
+12-month rolling mean, so monthly averaging discards very little.
+
+Months backed by fewer than half their days are dropped rather than reported.
+The record's first and last months are partial by construction, and the fixed
+365-day year rolls the final record into a January of its own — which would
+otherwise surface as a confident-looking monthly value resting on a single day.
+}
+
+\keyword{internal}
diff --git a/tests/testthat/test-climate-indices.R b/tests/testthat/test-climate-indices.R
index 690d740..8688280 100644
--- a/tests/testthat/test-climate-indices.R
+++ b/tests/testthat/test-climate-indices.R
@@ -19,6 +19,70 @@ psl_lines <- c(
   "  Calculated at NOAA PSL1"
 )
 
+# The retroflection index, in the shape Nature publishes it: two caption lines, a
+# header, then daily rows dated as decimal years on a fixed 365-day year.
+decimal_year_lines <- function(start = 1993, months = 1:12, days_in_last = 31) {
+  rows <- unlist(lapply(months, function(m) {
+    first <- as.Date(paste(start, m, 1, sep = "-"))
+    n <- if (m == max(months)) days_in_last else
+      as.integer(lubridate::days_in_month(first))
+    doy <- as.integer(format(first, "%j")) + seq_len(n) - 1
+    sprintf("%.11f,%.4f", start + (doy - 1) / 365, m / 100)
+  }))
+  c("Figure 3,", ",", "Date,Retroflection index", rows)
+}
+
+test_that("a decimal-year CSV parses to monthly means", {
+  series <- parse_index_table(decimal_year_lines(), format = "decimal_year_csv",
+                              name = "LCR")
+
+  expect_setequal(names(series), c("YEAR", "MONTH", "LCR"))
+  expect_equal(nrow(series), 12)
+  # Each month's daily values are all month/100, so the monthly mean is too.
+  expect_equal(series$LCR, (1:12) / 100)
+})
+
+test_that("a month backed by too few days is dropped, not reported", {
+  # December carries a single day here. Averaging it would produce a confident
+  # December value resting on one observation, which is the failure mode
+  # min_coverage guards against elsewhere in the package.
+  series <- parse_index_table(decimal_year_lines(days_in_last = 1),
+                              format = "decimal_year_csv", name = "LCR")
+
+  expect_equal(nrow(series), 11)
+  expect_false(12 %in% series$MONTH)
+
+  # Half a month is enough to keep.
+  half <- parse_index_table(decimal_year_lines(days_in_last = 16),
+                            format = "decimal_year_csv", name = "LCR")
+  expect_true(12 %in% half$MONTH)
+})
+
+test_that("decimal years are read on the 365-day convention the file uses", {
+  # The published file steps by exactly 1/365 between days. Reading it as 365.25
+  # would walk the derived dates off by several days across the record, moving
+  # values into neighbouring months.
+  series <- parse_index_table(decimal_year_lines(), format = "decimal_year_csv",
+                              name = "LCR")
+
+  expect_setequal(series$MONTH, 1:12)
+  expect_true(all(series$YEAR == 1993))
+})
+
+test_that("the retroflection index carries its citation", {
+  # It is the published output of one study rather than an operational product,
+  # so the reference has to travel with it.
+  entry <- climate_indices()$LCR
+
+  expect_match(entry$reference, "Jutras")
+  expect_match(entry$reference, "10\\.1038/s41467-023-38321-y")
+
+  dictionary <- as.data.frame(index_dictionary())
+  expect_match(dictionary$reference[dictionary$name == "LCR"], "Nature Communications")
+  # The operational indices have no single paper, and must not have one invented.
+  expect_true(all(is.na(dictionary$reference[dictionary$name %in% c("NAO", "AO")])))
+})
+
 test_that("a CPC table parses to one row per month", {
   series <- parse_index_table(cpc_lines, format = "cpc_table", name = "NAO")