From 8ae9b9281885850f6efb72cde2ac4900578ad7e0 Mon Sep 17 00:00:00 2001
From: sermomon <59225676+sermomon@users.noreply.github.com>
Date: Fri, 28 Apr 2023 09:35:24 +0200
Subject: [PATCH] Update kndvi.r

---
 kndvi.r | 54 +++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 45 insertions(+), 9 deletions(-)

diff --git a/kndvi.r b/kndvi.r
index df259d7..f1dd7d3 100644
--- a/kndvi.r
+++ b/kndvi.r
@@ -1,11 +1,47 @@
-# kNDVI with RBF kernel.
-# nir: near infrared (pixel) value
-# red: red value
+#if (!require("raster")) install.packages("raster")
+require(raster)
 
-# Fix or estimate a reasonable sigma value
-sigma <- 0.15
-knr <- exp(-(nir-red)^2/(2*sigma^2))
-kndvi <- (1-knr) / (1+knr)
+# Function to convert Sentinel-2 L2A into reflectance values
+s2_surface_reflectance <- function(x){
+  #@x: RasterLayer, RasterStack or RasterBrick containing Sentinel-2 L2A data
+  return(x/10000)
+}
 
-# An equivalent way to compute it is
-# kndvi <- tanh(((nir-red)/(2*sigma))^2)
+# Function to compute kNDVI from two RasterLayers
+kndvi <- function(nir, red, sigma=0.15, method="rbf", faster=FALSE){
+  #@nir: RasterLayer containing reflectance values in the Near-Infrared band
+  #@red: RasterLayer containing reflectance values in the Red band
+  #@sigma: sigma value. This argument accept a single sigma value or a RasterLayer containing the pixel-wise
+  # sigma value. By default sigma=0.15.
+  #@mode: method used to compute kNDVI. Options: "rbf" or "tanh". Default="rbf"
+  #@faster: Not implemented! If faster=TRUE non-evaluate range of values and make the function faster
+  if(max(raster::values(nir))>1){message("The NIR band contains values greater than 1")}
+  if(min(raster::values(nir))<0){message("The NIR band contains values lower than 0")}
+  if(max(raster::values(red))>1){message("The Red band contains values greater than 1")}
+  if(min(raster::values(red))<1){message("The Red band contains values lower than 0")}
+  if(method=="rbf"){
+    knr <- exp(-(nir-red)^2/(2*sigma^2))
+    kndvi <- (1-knr) / (1+knr)
+  }
+  if(method=="tanh"){
+    kndvi <- tanh(((nir-red)/(2*sigma))^2)
+  }
+  return(kndvi)
+}
+
+# Function to compute the pixel-wise sigma value
+pixel_wise_sigma <- function(nir, red){
+  #@nir: RasterLayer containing reflectance values in the Near-Infrared band
+  #@red: RasterLayer containing reflectance values in the Red band
+  sigma <- 0.5 * (nir + red)
+  return(sigma)
+}
+
+# Function to compute the median pixel-wise sigma value of the scene
+median_dist_sigma <- function(pixel_wise_sigma, nsamples=1000){
+  #@nir: RasterLayer containing reflectance values in the Near-Infrared band
+  #@red: RasterLayer containing reflectance values in the Red band
+  #@nsamples: Not implemented! Number of random samples used to compute the median value
+  set.seed(1234)
+  return(median(raster::values(pixel_wise_sigma)))
+}