Import Libraries

library(CubeR)
library(tibble)
library(magrittr)
library(raster)
library(tidyr)
library(lubridate)
library(purrr)
library(MonalisR)
library(dplyr)
library(ggplot2)

The norm_diff_pixel() function allows to calculate any Normalized Difference (e.g. NDVI or NDSI) between two bands. The following example extracts the NDVI from one Pixel over the complete Temporal Extent of the Data Cube.

coverage<-getCapability()[9]
bands <- coverage_get_bands(coverage = coverage)
coords<-c("620812.02", "5171499.36") %>% as.numeric
NIR <- bands[5]
Red <- bands[4]

nd1.tab <- norm_diff_pixel(coverage, coords, band1=NIR, band2=Red, date = NULL,plot = F) %>% as.tibble
#nd1.plot <- norm_diff_pixel(coverage, coords, band1=NIR, band2=Red, date = NULL,plot = T)

nd1.tab
## # A tibble: 171 x 2
##    times                        res
##    <fct>                      <dbl>
##  1 2015-06-27T00:00:00.000Z 0.498  
##  2 2015-07-04T00:00:00.000Z 0.459  
##  3 2015-07-24T00:00:00.000Z 0.799  
##  4 2015-08-03T00:00:00.000Z 0.924  
##  5 2015-08-06T00:00:00.000Z 0.908  
##  6 2015-08-13T00:00:00.000Z 0.905  
##  7 2015-08-16T00:00:00.000Z 0.0717 
##  8 2015-08-23T00:00:00.000Z 0.00936
##  9 2015-08-26T00:00:00.000Z 0.905  
## 10 2015-09-02T00:00:00.000Z 0.0261 
## # ... with 161 more rows
#plot(nd1.plot)

The norm_diff_raster() function allows to calculate any Normalized Difference (e.g. NDVI or NDSI) between two bands. It is similar to the previous function but works spatially at one date. It is possible to extract the whole image or only a subset

coverage<-getCapability()[9]
NIR <- bands[5]
Red <- bands[4]
bands <- coverage_get_bands(coverage = coverage)
coords<-c("620812.02", "5171499.36") %>% as.numeric
slice_E <- as.character(c(coords[1]-5000, coords[1]+5000))
slice_N <- as.character(c(coords[2]-5000, coords[2]+5000))

times <- coverage_get_timestamps(coverage = coverage)
date<-times[8]

nd1.tiff <- norm_diff_raster(coverage, slice_E, slice_N, date=date, band1=NIR, band2=Red)

nd1.tiff
## class       : RasterLayer 
## dimensions  : 1001, 1001, 1002001  (nrow, ncol, ncell)
## resolution  : 9.99001, 9.99001  (x, y)
## extent      : 615812, 625812, 5166499, 5176499  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32632 +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.4352941, 0.8980392  (min, max)

In the last examples we combined the functionality of NDVI computation provided by this package with an in-situ database implemented at EURAC. This so called MONALISA Database holds NDVI measurements in 15-minute temporal resolution and can be accessed with a second R-package developed at EURAC Research.

# Download Sentinel Data
coverage<-getCapability()[9] # S2 10m Reflectance
bands <- coverage_get_bands(coverage = coverage)
coords<-c("620812.02", "5171499.36") %>% as.numeric

nd1.tab <- norm_diff_pixel(coverage, coords, band1= bands[5], band2= bands[4], plot = F) %>% as.tibble

# Download InSitu Data
mnls_down<-downloadMonalisa(foi="vimes1500",
                            datestart = "2016-03-01 00:00",
                            dateend = "2016-12-31 00:00",
                            property = "Normalized Difference Vegetation Index - average",
                            procedure = "SpectralReflectanceSensor_vimes1500") %>% .[[1]] %>% as.tibble

# Format the Datasets
ndvi.MONA<-mnls_down %>%
  tidyr::separate(TimeStamp,c("Date","Time"),sep=" ") %>% 
  dplyr::filter(.,Time=="10:30:00") %>% 
  dplyr::select(-Time) %>% 
  dplyr::mutate(Date=as_date(Date)) %>% 
  tibble::add_column(Sensor="MONALISA NDVI")

ndvi.S2<-nd1.tab %>% 
  dplyr::mutate(Date=as_date(times)) %>% 
  tibble::add_column(Value=nd1.tab$res) %>% 
  tibble::add_column(Sensor="Raw S2 NDVI 10m") %>% 
  dplyr::select(-c(times,res)) %>% 
  dplyr::filter(year(Date)=="2016" & month(Date)>2)

# Bind
allndvi<-rbind(ndvi.MONA,ndvi.S2) %>% arrange(Date)
allndvi
## # A tibble: 364 x 3
##    Date         Value Sensor         
##    <date>       <dbl> <chr>          
##  1 2016-03-01 -0.0380 MONALISA NDVI  
##  2 2016-03-02  0.307  MONALISA NDVI  
##  3 2016-03-03 -0.0690 MONALISA NDVI  
##  4 2016-03-03  0.0186 Raw S2 NDVI 10m
##  5 2016-03-04 -0.0380 MONALISA NDVI  
##  6 2016-03-05 -0.0740 MONALISA NDVI  
##  7 2016-03-06 -0.0990 MONALISA NDVI  
##  8 2016-03-07 -0.0760 MONALISA NDVI  
##  9 2016-03-08 -0.0360 MONALISA NDVI  
## 10 2016-03-09 -0.0190 MONALISA NDVI  
## # ... with 354 more rows
ggplot(allndvi,aes(x=Date,y=Value,color=Sensor))+
  geom_line()+
  ggtitle("Combination of Ground and Satellite Data")+
  ylab("NDVI")