diff --git a/NDVI/README.md b/NDVI/README.md index 318dc702..e869ac8c 100644 --- a/NDVI/README.md +++ b/NDVI/README.md @@ -4,11 +4,11 @@ The data file (ndvi.csv) is composed of remotely sensed raw values of NDVI for a -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): -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: @@ -16,9 +16,37 @@ Other sources in the dataset include: * 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. + + + + + +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 + +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. diff --git a/NDVI/figures/ndvi-sensor-comparion-histograms.png b/NDVI/figures/ndvi-sensor-comparion-histograms.png new file mode 100644 index 00000000..35d3496f Binary files /dev/null and b/NDVI/figures/ndvi-sensor-comparion-histograms.png differ diff --git a/NDVI/figures/ndvi-sensor-comparion-one-to-one.png b/NDVI/figures/ndvi-sensor-comparion-one-to-one.png new file mode 100644 index 00000000..2e262c96 Binary files /dev/null and b/NDVI/figures/ndvi-sensor-comparion-one-to-one.png differ diff --git a/NDVI/portal-ndvi-shift.R b/NDVI/portal-ndvi-shift.R new file mode 100644 index 00000000..d20fbbf6 --- /dev/null +++ b/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)