Add information for correcting NDVI across sensors [no version bump]

NDVI/README.md

@@ -4,21 +4,49 @@ The data file (ndvi.csv) is composed of remotely sensed raw values of NDVI for a
 
 <img src="figures/portal_area.png" width="400px">
 
-The current data sources are Landsat 8 and 9. e.g. The Portal area is covered by 2 tiles each from the Landsat8 and 9 paths:
+The current data source is Landsat 9 and the Portal area is covered by 2 tiles from the Landsat path (Landsat 8 tiles also shown for comparison):
 
 <img src="figures/tiles.png" width="400px">
 
-The current source is the USGS/EROS api. NDVI is calculated as (B5 - B4)/(B5 + B4). Thus, the USGS source indicates a distinct time series from the AWS source, in which sr_ndvi was provided pre-calculated.
+The current source is the USGS/EROS API. NDVI is calculated as (B5 - B4)/(B5 + B4). Thus, the USGS source indicates a distinct time series from the AWS source, in which sr_ndvi was provided pre-calculated.
 
 Other sources in the dataset include:
 
 * Landsat (sensors 5-9, from a variety of sources)
 * GIMMS (an ensemble product from various AVHRR instruments on NOAA satellites)
 * MODIS (one instrument aboard the Terra & Aqua satellites)
 
-Further details on the efforts to create an ndvi time series at the site can be found at https://github.com/weecology/NDVIning.
+Due to the length of the the time-series we focus on the use of landsat data for analysis and a combined time-series across Landsats 5-9 is what is provided by default by the [`ndvi()`](https://weecology.github.io/portalr/reference/ndvi.html) function in our [portalr](https://weecology.github.io/portalr/) R package.
 
-### Other files
+## Correcting for differences between sensors
+
+The sensors on the different Landsat missions differ and this can produce shifts in the resulting NDVI values across the time-series.
+Each Landsat mission overlaps with the previous mission for the purpose of comparison (see [portal-ndvi-shift.R](portal-ndvi-shift.R)).
+We have compared each mission to the previous mission using data from the two sensors that is collected withing one of each other to evaluate these differences.
+
+<img src="figures/ndvi-sensor-comparison-one-to-one.png" width="400px">
+
+<img src="figures/ndvi-sensor-comparison-histograms.png" width="400px">
+
+Since notable shifts are present and they seem to be consistent across NDVI values we recommend correcting the means to a common mean based on Landsat 5.
+Calulated corrections to the previous sensor for Landsats 7-9 are:
+
+```
+# A tibble: 3 × 2
+  sensor2  correction_to_prev_sensor
+  <chr>                        <dbl>
+1 Landsat7                   0.00179
+2 Landsat8                  -0.0361
+3 Landsat9                  -0.00249
+```
+
+So to calculate a value for Landsat 9 one would first correct it to Landsat 8 and then to Landsat 7 and then to Landsat 5:
+
+> Landsat9_corrected = Landsat9_raw - 0.00249 - 0.361 + 0.00179
+
+Additional details of the history of satelittes and sensors along with earlier efforts to create a continuous NDVI time series for the site can be found at https://github.com/weecology/NDVIning.
+
+## Other files
 
 An older version of the ndvi time series (NDVI_monthly.csv) is retained for compatibility. It was obtained from the GIMMS NDVI3g dataset, which is a compilation of AVHRR weather sattelite data (Pinzon & Tucker 2014). The original data has a spatial resolution of 8km and a twice monthly temporal resolution, with availability from 1981-2013. It is averaged to monthly NDVI values. 
 

NDVI/portal-ndvi-shift.R

@@ -0,0 +1,89 @@
+library(dplyr)
+library(ggplot2)
+library(lubridate)
+library(portalr)
+library(purrr)
+library(readr)
+library(stringr)
+
+find_sensor_matches <- function(ndvi_data, max_days = 10) {
+  ndvi_data$date <- as.Date(ndvi_data$date)
+  sensors <- c("Landsat5", "Landsat7", "Landsat8", "Landsat9")
+  results_list <- list()
+
+  for (sensor_i in 1:(length(sensors) - 1)) {
+    sensor1_data <- ndvi_data |> filter(sensor == sensors[sensor_i])
+    sensor2_data <- ndvi_data |> filter(sensor == sensors[sensor_i + 1])
+    for (i in 1:nrow(sensor1_data)) {
+      current_date <- sensor1_data$date[i]
+      current_ndvi <- sensor1_data$ndvi[i]
+      current_pixel_count <- sensor1_data$pixel_count[i]
+      date_diffs <- abs(as.numeric(sensor2_data$date - current_date))
+      min_diff_idx <- which.min(date_diffs)
+      min_diff_days <- date_diffs[min_diff_idx]
+      if (min_diff_days <= max_days) {
+        match_row <- sensor2_data[min_diff_idx, ]
+        result_row <- tibble(
+          sensor1 = sensors[sensor_i],
+          ndvi1 = current_ndvi,
+          date1 = current_date,
+          pixel_count1 = current_pixel_count,
+          sensor2 = sensors[sensor_i + 1],
+          ndvi2 = match_row$ndvi,
+          date2 = match_row$date,
+          pixel_count2 = match_row$pixel_count,
+          days_diff = min_diff_days
+        )
+        results_list <- append(results_list, list(result_row))
+      }
+    }
+  }
+  results <- bind_rows(results_list)
+  return(results)
+}
+
+ndvi <- load_datafile("NDVI/ndvi.csv") |>
+  as_tibble()
+
+sensor_matches <- find_sensor_matches(ndvi, max_days = 10)
+
+ndvi_pairs <- sensor_matches |>
+  filter(
+    days_diff <= 10,
+    ndvi1 > 0,
+    ndvi2 > 0,
+    pixel_count1 > 3562 / 2,
+    pixel_count2 > 3562 / 2,
+  )
+
+date_ranges <- ndvi |>
+  filter(str_detect(sensor, "Landsat")) |>
+  group_by(sensor) |>
+  summarise(
+    min_date = min(date, na.rm = TRUE),
+    max_date = max(date, na.rm = TRUE),
+    .groups = "drop"
+  )
+
+print(date_ranges)
+
+ggplot(ndvi_pairs, aes(x = ndvi1, y = ndvi2)) +
+  geom_point() +
+  geom_abline(slope = 1, intercept = 0) +
+  facet_wrap(vars(sensor1, sensor2), scales = "free")
+
+ggsave("ndvi-sensor-comparion-one-to-one.png")
+
+ggplot(ndvi_pairs, (aes(x = ndvi1 - ndvi2))) +
+  geom_histogram() +
+  geom_vline(xintercept = 0) +
+  xlab("Sensor Difference (Current - Previous)") +
+  facet_wrap(vars(sensor2), ncol = 1)
+
+ggsave("ndvi-sensor-comparion-histograms.png")
+
+instrument_diffs <- ndvi_pairs |>
+  group_by(sensor2) |>
+  summarize(correction_to_prev_sensor = mean(ndvi1 - ndvi2))
+
+print(instrument_diffs)