From 72cd7cce72e1454566bae2a551ebbf73bef1507f Mon Sep 17 00:00:00 2001 From: David Klehr Date: Wed, 14 Aug 2024 09:39:48 +0200 Subject: [PATCH 1/4] Added Spline Reconstruction --- .../Bolten_Spline_FORCE_UDF.R | 135 ++++++++++++++++++ rstats/ts/spline-reconstruction/readme.md | 20 +++ 2 files changed, 155 insertions(+) create mode 100644 rstats/ts/spline-reconstruction/Bolten_Spline_FORCE_UDF.R create mode 100644 rstats/ts/spline-reconstruction/readme.md diff --git a/rstats/ts/spline-reconstruction/Bolten_Spline_FORCE_UDF.R b/rstats/ts/spline-reconstruction/Bolten_Spline_FORCE_UDF.R new file mode 100644 index 0000000..d40e3a4 --- /dev/null +++ b/rstats/ts/spline-reconstruction/Bolten_Spline_FORCE_UDF.R @@ -0,0 +1,135 @@ +force_rstats_init <- function(dates, sensors, bandnames){ + + # Year which should be reconstructed + year_to_interpolate <- 2023 + # Days to predict in this year and the intervall: 60 to 330 (1st March to 26th November) + DOYs_to_predict <- seq(60,330,by =10) + dates_to_predict <- as.Date(paste(year_to_interpolate, DOYs_to_predict), format = "%Y %j") + + band_names <- paste(substr(as.character(dates_to_predict),1,4),substr(as.character(dates_to_predict),6,7),substr(as.character(dates_to_predict),9,10), sep ="") + return(band_names) +} + +force_rstats_pixel <- function(inarray, dates, sensors, bandnames, nproc){ + + # Year which should be interpolated (same like above) + year_to_interpolate <- 2023 + # Days to predict in this yearand the intervall: 60 to 330 (1st March to 26th November) + DOYs_to_predict <- seq(60,330,by =10) + dates_to_predict <- as.Date(paste(year_to_interpolate, DOYs_to_predict), format = "%Y %j") + + # spline variables + # smoothing factor for the spline reconstruction + smooth_fac <- 0.5 + # Bolton's variable of maximum weight to assing for the predessesor years + # the year of reconstruction has a wheight of 1 + max_weight <- 0.2 + + # cumulate the DOY to the year of interpolation + # start year 2015 (example), because of Sentinel 2 launch date, for e.g. Landsat adjust to your time span + start_year <- 2015 + DOYs_to_predict <- (year_to_interpolate - start_year) * 365 + DOYs_to_predict + + tryCatch({ + # grap FORCE no-data case + if (all(is.na(inarray[,1]))){ + return(rep(-9999,length(DOYs_to_predict))) + } + + # calculate cumulative DOYs for the input data + DOYs <- as.numeric(format(dates, "%j")) + years <- as.numeric(substr(as.character(dates),1,4)) + cumulative_DOYs <- (years - start_year) * 365 + DOYs + + # join the data to a dataframe + df <- data.frame(x=cumulative_DOYs,y=inarray[,1]) + + # ------- 1.1 calcualte Mean Function -------------- + euc.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2)) + + # define Start and endpoints for the three spline reconstuctions + DOY_borders_year <- c((year_to_interpolate-start_year)*365 - 180, (year_to_interpolate-start_year+1)*365+180) + DOY_borders_b <- c((year_to_interpolate-start_year-1)*365 - 180, (year_to_interpolate-start_year)*365+180) + DOY_borders_bb <- c((year_to_interpolate-start_year-2)*365 - 180, (year_to_interpolate-start_year-1)*365+180) + + # create dataframes for reconstuction + data_year <- na.exclude(df[df$x %in% seq(DOY_borders_year[1],DOY_borders_year[2]),c(1,2)]) + data_b <- na.exclude(df[df$x %in% seq(DOY_borders_b[1] ,DOY_borders_b[2]),c(1,2)]) + data_bb <- na.exclude(df[df$x %in% seq(DOY_borders_bb[1],DOY_borders_bb[2]),c(1,2)]) + + # calculate spline model for year of reconstruction and predict + DOYs_target_year <- seq(DOY_borders_year[1],DOY_borders_year[2]) + tryCatch({ + eval( spline_model_year <<- smooth.spline(data_year$x,data_year$y, spar =smooth_fac) ) + eval( predict_year <<- predict(spline_model_year, x = DOYs_target_year)$y ) + }, error = function(err) { + return(rep(-9999,length(DOYs_to_predict))) + }) + + #calculate d_max + mean_year <- mean(na.exclude(data_year$y)) + mean_prediction <- rep(mean_year,length(DOYs_target_year)) + d_max = euc.dist(mean_prediction, predict_year) / 10000 + + + # --------- 1.2 spline for precessor years ------------ + # one year before + # predict with DOYs of year of reconstruction, for difference calculation + # between the two spline reconstructions + tryCatch({ + eval( spline_model_b <<- smooth.spline(data_b$x+365,data_b$y, spar =smooth_fac) ) + eval( predict_b <<- predict(spline_model_b, x = DOYs_target_year)$y ) + }, error = function(err) { + return(rep(-9999,length(DOYs_to_predict))) + }) + d_b = euc.dist(predict_year, predict_b) / 10000 + + # two years before + # predict with DOYs of year of reconstruction, for difference calculation + # between the two spline reconstructions + tryCatch({ + eval( spline_model_bb <<- smooth.spline(data_bb$x+(365*2),data_bb$y, spar =smooth_fac) ) + eval( predict_bb <<- predict(spline_model_bb, x = DOYs_target_year)$y ) + }, error = function(err) { + return(rep(-9999,length(DOYs_to_predict))) + }) + d_bb = euc.dist(predict_year, predict_bb) / 10000 + + # ---------- 1.3 Calculate weights ------------------- + # one year before + if (d_max != 0) { + weight_b = max_weight*(1-(d_b/d_max)) + if (weight_b < 0){ + weight_b = 0 + } + } else {weight_b = 0} + + # two years before + if (d_max != 0) { + weight_bb = max_weight*(1-(d_bb/d_max)) + if (weight_bb < 0){ + weight_bb = 0 + } + } else {weight_bb = 0} + + #----------- 1.4 apply weights and calculate weighted spline -------------- + # combine the time series to one year and assign weights to the new data points + combined_x <- c(data_year$x , (data_b$x+365)[weight_b>0] , (data_bb$x + (365*2))[weight_bb>0]) + combined_y <- c(data_year$y , data_b$y[weight_b>0] , data_bb$y[weight_bb>0]) + vec_weights <- c(rep(1,length(data_year$x)), + rep(weight_b,length(data_b$x))[weight_b>0], + rep(weight_bb,length(data_bb$x))[weight_bb>0]) + + # calculate weighted spline + tryCatch({ + eval( spline_model_combined <<- smooth.spline(x=combined_x, y=combined_y, w = vec_weights , spar =smooth_fac) ) + eval( predict_combined <<- predict(spline_model_combined, x = DOYs_to_predict)$y ) + }, error = function(err) { + return(rep(-9999,length(DOYs_to_predict))) + }) + return(predict_combined) + + }, error = function(err) { + return(rep(-9999,length(DOYs_to_predict))) + }) +} \ No newline at end of file diff --git a/rstats/ts/spline-reconstruction/readme.md b/rstats/ts/spline-reconstruction/readme.md new file mode 100644 index 0000000..5b96f6a --- /dev/null +++ b/rstats/ts/spline-reconstruction/readme.md @@ -0,0 +1,20 @@ +# Time series reconstruction for forest using splines + +© +Copyright 2024, David Klehr + +## Run with + +- program: ``force-higher-level`` +- submodule: ``TSA`` +- DATE_RANGE: ``xxxx-07-01 yyyy-06-31`` +xxxx = three years before your target year +yyyy = one year after your target year +e.g. for target year 2022: 2019-07-01 2023-06-31 +- UDF type: ``RSTATS_TYPE = PIXEL`` +- required parameters:``none`` +- required R libraries: ``numpy, scipy`` + +## References + +- Bolton, D.K., Gray, J.M., Melaas, E.K., Moon, M., Eklundh, L., Friedl, M.A., 2020. **Continental-scale land surface phenology from harmonized Landsat 8 and Sentinel-2 imagery**. *Remote Sensing of Environment 240*, 111685. https://doi.org/10.1016/j.rse.2020.111685. \ No newline at end of file From eef67b8c0f35322c89232ea8f538e4d2096529a2 Mon Sep 17 00:00:00 2001 From: David Klehr Date: Wed, 14 Aug 2024 09:41:16 +0200 Subject: [PATCH 2/4] changes in readme --- rstats/ts/spline-reconstruction/readme.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rstats/ts/spline-reconstruction/readme.md b/rstats/ts/spline-reconstruction/readme.md index 5b96f6a..e03de31 100644 --- a/rstats/ts/spline-reconstruction/readme.md +++ b/rstats/ts/spline-reconstruction/readme.md @@ -8,12 +8,13 @@ Copyright 2024, David Klehr - program: ``force-higher-level`` - submodule: ``TSA`` - DATE_RANGE: ``xxxx-07-01 yyyy-06-31`` + xxxx = three years before your target year yyyy = one year after your target year e.g. for target year 2022: 2019-07-01 2023-06-31 - UDF type: ``RSTATS_TYPE = PIXEL`` - required parameters:``none`` -- required R libraries: ``numpy, scipy`` +- required R libraries: ``none`` ## References From 2d58980373271e0aa2215751a7f10dca6d7cf004 Mon Sep 17 00:00:00 2001 From: David Klehr Date: Wed, 14 Aug 2024 09:43:13 +0200 Subject: [PATCH 3/4] changes readme --- rstats/ts/spline-reconstruction/readme.md | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/rstats/ts/spline-reconstruction/readme.md b/rstats/ts/spline-reconstruction/readme.md index e03de31..cbb83d2 100644 --- a/rstats/ts/spline-reconstruction/readme.md +++ b/rstats/ts/spline-reconstruction/readme.md @@ -8,10 +8,9 @@ Copyright 2024, David Klehr - program: ``force-higher-level`` - submodule: ``TSA`` - DATE_RANGE: ``xxxx-07-01 yyyy-06-31`` - -xxxx = three years before your target year -yyyy = one year after your target year -e.g. for target year 2022: 2019-07-01 2023-06-31 + * xxxx = three years before your target year + * yyyy = one year after your target year + * e.g. for target year 2022: 2019-07-01 2023-06-31 - UDF type: ``RSTATS_TYPE = PIXEL`` - required parameters:``none`` - required R libraries: ``none`` From ea53aebe5a080b2e8721e24a9d4ffd04f474a957 Mon Sep 17 00:00:00 2001 From: David Klehr Date: Wed, 14 Aug 2024 09:43:51 +0200 Subject: [PATCH 4/4] changes readme --- rstats/ts/spline-reconstruction/readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rstats/ts/spline-reconstruction/readme.md b/rstats/ts/spline-reconstruction/readme.md index cbb83d2..c03e0ff 100644 --- a/rstats/ts/spline-reconstruction/readme.md +++ b/rstats/ts/spline-reconstruction/readme.md @@ -10,7 +10,7 @@ Copyright 2024, David Klehr - DATE_RANGE: ``xxxx-07-01 yyyy-06-31`` * xxxx = three years before your target year * yyyy = one year after your target year - * e.g. for target year 2022: 2019-07-01 2023-06-31 + * e.g. for target year 2022: ``2019-07-01 2023-06-31`` - UDF type: ``RSTATS_TYPE = PIXEL`` - required parameters:``none`` - required R libraries: ``none``