| Title: | Vegetation Phenological Cycle and Anomaly Detection using Remote Sensing Data |
|---|---|
| Description: | Calculates phenological cycle and anomalies using a non-parametric approach applied to time series of vegetation indices derived from remote sensing data or field measurements. The package implements basic and high-level functions for manipulating vector data (numerical series) and raster data (satellite derived products). Processing of very large raster files is supported. For more information, please check the following paper: Chávez et al. (2023) <doi:10.3390/rs15010073>. |
| Authors: | Roberto O. Chávez [cre, aut], Sergio A. Estay [cre, aut], José A. Lastra [cre,ctb] & Carlos G. Riquelme [ctb] |
| Maintainer: | José A. Lastra <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 2.0.0 |
| Built: | 2025-01-27 06:11:42 UTC |
| Source: | https://github.com/labgrs/npphen |
Based on the annual reference frequency distribution (rfd) of a vegetation index time series (e.g. a numeric vector of NDVI), calculates anomalies and how extreme these anomalies are (rfd position ranging from 0 to 100)
ExtremeAnom(x, dates, h, refp, anop, rge, output = "both", rfd = 0.9)ExtremeAnom(x, dates, h, refp, anop, rge, output = "both", rfd = 0.9)
x |
Numeric vector. A time series of a vegetation index (e.g. LAI, NDVI, EVI) or any other variable with seasonal behavior. The code has been optimized to work with integer values. Please re-scale the input values if necessary (e.g. NDVI ranging from 0.0000 to 1.0000, multiply by 10,000 |
dates |
A date vector. The number of dates must be equal to the number of "x" values (numeric input vector). |
h |
Numeric. Indicates the geographic hemisphere. This argument defines the starting date of the growing season. h = 1 for the Northern Hemisphere (season starting on January 1st), h = 2 for the Southern Hemisphere (season starting on July 1st). |
refp |
Numeric vector with the correlative number of dates to be used as reference period. For example, refp = c(1:388) for MODIS Vegetation Index 16-days composites MOD13Q1 (18/02/2000 – 18/12/2016) |
anop |
Numeric vector with the correlative number of dates for the period in which the anomalies and rfd position (how extreme the anomalies are) will be calculated. For example anop = c(389:411) for the complete 2017 of MODIS Vegetation Index 16-days composites MOD13Q1 (01/01/2017 – 19/12/2017). anop y refp can be overlapped. |
rge |
Numeric vector with two values setting the minimum and maximum values of the response variable (e.g. NDVI) used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge = c(0,10000) |
output |
Character string. Defines the output values. 'both' (default) returns both anomalies and rfd position together as a single numeric vector, 'anomalies' returns only anomalies, 'rfd' returns only rfd values (how extreme the anomalies are) and 'clean' returns only extreme anomalies, i.e. anomalies at which a given rfd is overpass (e.g. 0.90). This critical threshold is set by the users using the rfd argument. |
rfd |
Numeric. This argument only applies when the argument output = 'clean'. It defines the percentile (from 0 to 0.99) of the reference frequency distribution, for which anomalies are not flagged as extreme anomalies. For example, if 'rfd = 0.90' only anomalies falling outside the '0.90 rfd' (default value) will be flagged as extreme anomalies while the rest will be neglected (NA values). Please notice that 'rfd = 0.90' implies that the 5% of the most extreme positive and 5% of the most extreme negative anomalies will be considered. |
Calculates anomalies and a probabilistic measure of how extreme the anomalies are using a numeric vector of vegetation canopy greenness, e.g. satellite based Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). For this purpose, it divides the time series (numeric vector) of vegetation greenness into 2: the reference period, from which the annual phenological reference is calculated (same as Phen function), and the observation period, for which we want to calculate anomalies. This anomalies can be filtered by the position of the observation within the historical rfd. Users can, for example, set 'rfd = 0.95' to consider only anomalies that outside the 95% rfd of historical records.
numerical vector
library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series using a raster stack and a date database from Central Chile # Obtain data from a particular pixel generating a time series md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Anomaly detection for the given pixel anomRFD_MD <- ExtremeAnom( x = md_pixelts, dates = modis_dates, h = 2, refp = c(1:423), anop = c(884:929), rfd = 0.9, output = "both", rge = c(0, 10000) ) # Basic plot selection <- which(anomRFD_MD[47:92] > 90) barplot( names = format.Date(modis_dates[884:929], format = "%Y-%m"), anomRFD_MD[1:46], col = ifelse(anomRFD_MD[1:46] < 0, "red", "blue"), main = "Anomalies whit rfd > 0.90" ) abline(h = 0) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), # showing extreme positive anomalies (greening)## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series using a raster stack and a date database from Northern Chile # Obtain data from a particular pixel generating a time series bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Anomaly detection for the given pixel anomRFD_BD <- ExtremeAnom( x = bd_pixelts, dates = modis_dates, h = 2, refp = c(1:423), anop = c(723:768), rfd = 0.9, output = "both", rge = c(0, 10000) ) # Basic plot selection <- which(anomRFD_BD[47:92] > 95) # basic plot barplot( names = format.Date(modis_dates[723:768], format = "%Y-%m"), anomRFD_BD[1:46], col = ifelse(anomRFD_BD[1:46] < 0, "red", "blue"), main = "Anomalies whit rfd > 0.95" ) abline(h = 0)library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series using a raster stack and a date database from Central Chile # Obtain data from a particular pixel generating a time series md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Anomaly detection for the given pixel anomRFD_MD <- ExtremeAnom( x = md_pixelts, dates = modis_dates, h = 2, refp = c(1:423), anop = c(884:929), rfd = 0.9, output = "both", rge = c(0, 10000) ) # Basic plot selection <- which(anomRFD_MD[47:92] > 90) barplot( names = format.Date(modis_dates[884:929], format = "%Y-%m"), anomRFD_MD[1:46], col = ifelse(anomRFD_MD[1:46] < 0, "red", "blue"), main = "Anomalies whit rfd > 0.90" ) abline(h = 0) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), # showing extreme positive anomalies (greening)## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series using a raster stack and a date database from Northern Chile # Obtain data from a particular pixel generating a time series bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Anomaly detection for the given pixel anomRFD_BD <- ExtremeAnom( x = bd_pixelts, dates = modis_dates, h = 2, refp = c(1:423), anop = c(723:768), rfd = 0.9, output = "both", rge = c(0, 10000) ) # Basic plot selection <- which(anomRFD_BD[47:92] > 95) # basic plot barplot( names = format.Date(modis_dates[723:768], format = "%Y-%m"), anomRFD_BD[1:46], col = ifelse(anomRFD_BD[1:46] < 0, "red", "blue"), main = "Anomalies whit rfd > 0.95" ) abline(h = 0)
Based on the annual reference frequency distribution (RFD) of a vegetation index time series (e.g. a raster stack of NDVI), calculates anomalies and how extreme these anomalies are (RFD position ranging from 0 to 100).
ExtremeAnoMap(s, dates, h, refp, anop, rge, output = "both", rfd = 0.9, nCluster, outname, datatype)ExtremeAnoMap(s, dates, h, refp, anop, rge, output = "both", rfd = 0.9, nCluster, outname, datatype)
s |
SpatRaster. A time series of a vegetation index (e.g. LAI, NDVI, EVI) or any other variable with seasonal behavior.The code has been optimized to work with integer values. Please re-scale the input SpatRaster if necessary (e.g. NDVI ranging from 0.0000 to 1.0000, multiply by 10,000). |
dates |
A date vector. The number of dates must be equal to the number of layers of the vegetation index SpatRaster. |
h |
Numeric. Indicates the geographic hemisphere to define the starting date of the growing season. h = 1 if the vegetation is in the Northern Hemisphere (season starting on January 1st), h = 2 if it is in the Southern Hemisphere (season starting on July 1st). |
refp |
Numeric vector with the correlative number of dates to be used as reference period. For example, refp = c(1:393) for MODIS Vegetation Index 16-days composites (18/02/2000 – 06/06/2017) |
anop |
Numeric vector with the correlative number of dates for the period in which the anomalies will be calculated. For example refp = c(21:43) for the first complete year for MODIS Vegetation Index 16-days composites (01/01/2001 – 19/12/2001). anop y refp can overlap |
rge |
Numeric vector with the minimum and maximum values of the vegetation index used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge = c(0,10000). |
output |
Character string. Defines the outputs. output = 'both' (default) returns both the vegetation index anomalies and rfd position together as a single numeric vector, output = 'anomalies' returns only anomalies, output = 'rfd' returns only rfd values (how extreme the anomalies are) and output = 'clean' returns only anomalies at which a given rfd is overpass (e.g. 0.90). This critical threshold is set by the users using the rfd argument. |
rfd |
Numeric. This argument only applies when the argument output = 'clean'. It defines the percentile (from 0 to 0.99) of the reference frequency distribution, for which anomalies are not flagged as extreme anomalies. For example, if 'rfd = 0.90' only anomalies falling outside the '0.90 rfd' (default value) will be flagged as extreme anomalies while the rest will be neglected (NA values). Please notice that'rfd = 0.90' implies that the 5% of the most extreme positive and 5% of the most extreme negative anomalies will be considered. |
nCluster |
Numeric. Number of CPU cores to be used for computational calculations |
outname |
Character vector with the output directory path and filename with extension or only the filename and extension if the working directory was set. For example outname = "output_anom.tif". See |
datatype |
Character. Output data type. For example datatype = "INT2S" for signed integer values. See |
Similar to ExtremeAnom, it calculates phenological anomalies but using a raster stack instead of a numeric vector of vegetation canopy greenness values (e.g. Leaf Area Index, LAI) or satellite based greenness proxies such as the Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). For this purpose, it divides the time series (raster stack) of vegetation greenness into 2: the reference period, from which the annual phenological cycle is calculated (same as PhenMap function), and the observation period, for which we want to calculate anomalies with respect to the annual phenological cycle. Negative anomalies correspond to observed values lower than the reference and positive anomalies to values higher than the reference. his anomalies can be filtered by the position of the observation within the historical rfd. Users can, for example, set 'rfd = 0.95' to consider only anomalies that outside the 95% rfd of historical records.
SpatRaster
## Not run: ## DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING## ## Testing with the MegaDrought_stack from Central Chile (NDVI), # showing extreme negative anomalies (browning)## # Load data # SpatRaster library(terra) f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Define the number of cores to be use. In this example we use 1 nc1 <- 1 ExtremeAnoMap( s = MegaDrought, dates = modis_dates, h = 2, refp = c(1:423), anop = c(884:929), rfd = 0.9, output = "both", nCluster = nc1, outname = "anomRFD_MD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map_an1 <- rast("anomRFD_MD.tif")#run only for load anomaly brick # plot(map_an1) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), # showing extreme positive anomalies (greening)## # Load data SpatRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Dates data("modis_dates") # Define the number of cores to be use. In this example we use 1 nc1 <- 1 ExtremeAnoMap( s = Bdesert, dates = modis_dates, h = 2, refp = c(1:423), anop = c(723:768), rfd = 0.9, output = "both", nCluster = nc1, outname = "anomRFD_BD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map_an1 <- rast("anomRFD_BD.tif")#run only for load anomaly brick # plot(map_an1) ## End(Not run)## Not run: ## DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING## ## Testing with the MegaDrought_stack from Central Chile (NDVI), # showing extreme negative anomalies (browning)## # Load data # SpatRaster library(terra) f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Define the number of cores to be use. In this example we use 1 nc1 <- 1 ExtremeAnoMap( s = MegaDrought, dates = modis_dates, h = 2, refp = c(1:423), anop = c(884:929), rfd = 0.9, output = "both", nCluster = nc1, outname = "anomRFD_MD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map_an1 <- rast("anomRFD_MD.tif")#run only for load anomaly brick # plot(map_an1) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), # showing extreme positive anomalies (greening)## # Load data SpatRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Dates data("modis_dates") # Define the number of cores to be use. In this example we use 1 nc1 <- 1 ExtremeAnoMap( s = Bdesert, dates = modis_dates, h = 2, refp = c(1:423), anop = c(723:768), rfd = 0.9, output = "both", nCluster = nc1, outname = "anomRFD_BD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map_an1 <- rast("anomRFD_BD.tif")#run only for load anomaly brick # plot(map_an1) ## End(Not run)
A vector with dates for both RasterStack datasets: Mega drought and Blooming Desert, with 929 dates
modis_datesmodis_dates
Vector of dates
The functions in this package estimate the expected annual phenological cycle from time series or raster stack of vegetation (greenness) indexes. The algorithm to estimate the annual phenological cycle (used by the functions in npphen) uses a bivariate kernel density estimator in the index-time space. In this space, the x-axis corresponds to days of the growing season (1-365) and the y-axis to the vegetation index values, which range's values are set using the rge argument (see each function's vignette for details). The expected value of the index for each day of the growing season (the expected phenology) is approximated by the maximum value of the kernel at that day. Anomalies are calculated as deviations from the expected values at given days. The package implements basic and high-level functions for manipulating vector data (numerical series) and raster data (satellite derived products). Processing of very large raster files is supported. For methodological details of kernel density estimation see Wand & Jones (1994).
Roberto O. Chávez [email protected]
Sergio A. Estay [email protected]
José A. Lastra [email protected]
Carlos G. Riquelme [email protected]
Wand, M.P. & Jones, M.C. (1994) Kernel smoothing. Crc Press.
Phen, PhenMap, ExtremeAnom, ExtremeAnoMap, PhenKplot
Estimates the annual phenological cycle from a time series of vegetation greenness.
Phen(x, dates, h, frequency = "16-days", rge, plot = T)Phen(x, dates, h, frequency = "16-days", rge, plot = T)
x |
Numeric vector. A time series of a vegetation index (e.g. LAI, NDVI, EVI) or any other variable with seasonal behavior. The code has been optimized to work with integer values. Please re-scale the input values if necessary (e.g. NDVI ranging from 0.0000 to 1.0000, multiply by 10,000 |
dates |
A date vector. The number of dates must be equal to the number of "x" values (numeric input vector). |
h |
Numeric. Indicates the geographic hemisphere to define the starting date of the growing season. h = 1 if the vegetation is in the Northern Hemisphere (season starting at January 1st), h = 2 if it is in the Southern Hemisphere (season starting at July 1st) |
frequency |
Character string. Defines the number of samples for the output phenology and must be one of the following options: 'daily' giving output vector of length 365, '8-days' giving output vector of length 46 (i.e MOD13Q1 and MYD13Q1), 'monthly' giving output vector of length 12,'bi-weekly' giving output vector of length 24 (i.e. GIMMS) or '16-days' (default) giving output vector of length 23 (i.e MOD13Q1 or MYD13Q1). |
rge |
Numeric vector with two values setting the minimum and maximum values of the response variable (e.g. NDVI) used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge = c(0,10000) |
plot |
Logical. Set TRUE (default) or FALSE to show the phenology curve plot. |
Derives the annual phenological cycle for a standard growing season using a numeric vector of vegetation canopy greenness values (e.g. Leaf Area Index, LAI) or satellite based greenness proxies such as the Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). A vector with dates for the greenness values is also required.
A numeric vector, where each value represents the expected greeness at that date
library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series using a raster stack and a date database from Central Chile # Obtain data from a particular pixel generating a time series md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Phenology for the given pixel Phen(x = md_pixelts, dates = modis_dates, h = 2, frequency = "16-days", rge = c(0, 10000)) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), h=2 ## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series using a raster stack and a date database from Northern Chile # Obtain data from a particular pixel generating a time series bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Phenology for the given pixel Phen(x = bd_pixelts, dates = modis_dates, h = 2, frequency = "16-days", rge = c(0, 10000))library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series using a raster stack and a date database from Central Chile # Obtain data from a particular pixel generating a time series md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Phenology for the given pixel Phen(x = md_pixelts, dates = modis_dates, h = 2, frequency = "16-days", rge = c(0, 10000)) ## Testing with the Bdesert_stack from the Atacama Desert, Northern Chile (NDVI), h=2 ## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series using a raster stack and a date database from Northern Chile # Obtain data from a particular pixel generating a time series bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Phenology for the given pixel Phen(x = bd_pixelts, dates = modis_dates, h = 2, frequency = "16-days", rge = c(0, 10000))
Plot the most probable vegetation greenness values.
PhenKplot(x, dates, h, xlab, ylab, rge)PhenKplot(x, dates, h, xlab, ylab, rge)
x |
Numeric vector. A time series of a vegetation index (e.g. LAI, NDVI, EVI) or any other variable with seasonal behavior. The code has been optimized to work with integer values. Please re-scale the input values if necessary (e.g. NDVI ranging from 0.0000 to 1.0000, multiply by 10,000) |
dates |
A date vector. The number of dates must be equal to the number of "x" values (numeric input vector). |
h |
Numeric. Indicates the geographic hemisphere to define the starting date of the growing season. h = 1 if the vegetation is in the Northern Hemisphere (season starting on January 1st), h = 2 if it is in the Southern Hemisphere (season starting on July 1st) |
xlab |
Character vector (or expression) giving plot title in x axis label (e.g. xlab = "day of the growing season") |
ylab |
Character vector (or expression) giving plot title in y axis label (e.g. ylab = "NDVI") |
rge |
Numeric vector with the minimum and maximum values of the vegetation index (e.g. NDVI) used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge = c(0,10000) |
It is a "heatmap" of the annual phenological variability of the Phen output. It calculates and plot a likelihood map of the vegetation-index–time space using a numeric vector of greenness proxies such as the Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). Also a vector with dates for the vegetation index values is required. This function is partially based on the ci2d function on package gplots.
library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series from a particular pixel # using a SpatRaster and a date for Central Chile md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Variability of the annual phenology for the given pixel PhenKplot(x = md_pixelts, dates = modis_dates, h = 2, xlab = "DGS", ylab = "NDVI", rge = c(0, 10000)) ## Testing with the Bdesert_spatRast from ## the Atacama Desert, Northern Chile (NDVI), h=2 ## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series from a particular pixel # using a SpatRaster and a date for Northern Chile bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Variability of the annual phenology for the given pixel PhenKplot(x = bd_pixelts, dates = modis_dates, h = 2, xlab = "DGS", ylab = "NDVI", rge = c(0, 10000))library(lubridate) library(terra) ## Testing raster data from Central Chile (NDVI), h=2## # Load data f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Generate a Raster time series from a particular pixel # using a SpatRaster and a date for Central Chile md_pixel <- cellFromXY(MegaDrought, cbind(313395, 6356610)) md_pixelts <- as.numeric(MegaDrought[md_pixel]) plot(modis_dates, md_pixelts, type = "l") # Variability of the annual phenology for the given pixel PhenKplot(x = md_pixelts, dates = modis_dates, h = 2, xlab = "DGS", ylab = "NDVI", rge = c(0, 10000)) ## Testing with the Bdesert_spatRast from ## the Atacama Desert, Northern Chile (NDVI), h=2 ## # Load data # SparRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Generate a Raster time series from a particular pixel # using a SpatRaster and a date for Northern Chile bd_pixel <- cellFromXY(Bdesert, cbind(286638, 6852107)) bd_pixelts <- as.numeric(Bdesert[bd_pixel]) plot(modis_dates, bd_pixelts, type = "l") # Variability of the annual phenology for the given pixel PhenKplot(x = bd_pixelts, dates = modis_dates, h = 2, xlab = "DGS", ylab = "NDVI", rge = c(0, 10000))
Estimates annual Land Surface Phenology (LSP) using time series of a vegetation greenness raster stack.
PhenMap(s, dates, h, frequency = "16-days", nCluster, outname, datatype, rge)PhenMap(s, dates, h, frequency = "16-days", nCluster, outname, datatype, rge)
s |
SpatRaster of a vegetation index (e.g. LAI, NDVI, EVI) or any other variable with seasonal behavior. The code has been optimized to work with integer values. Please re-scale the input SpatRaster if necessary (e.g. NDVI ranging from 0.0000 to 1.0000, multiply by 10,000) |
dates |
A date vector. The number of dates must be equal to the number of layers of the vegetation index SpatRaster. |
h |
Numeric. Indicates the geographic hemisphere to define the starting date of the growing season. h = 1 if the vegetation is in the Northern Hemisphere (season starting on January 1st), h = 2 if it is in the Southern Hemisphere (season starting on July 1st) |
frequency |
Character string. Defines the number of samples for the output phenology and must be one of the following options: 'daily' with an output vector of length 365, '8-days' with an output vector of length 46 (e.g. MOD13Q1 and MYD13Q1 combined), 'monthly' with an output vector of length 12,'bi-weekly' with an output vector of length 24 (i.e. GIMMS) or '16-days' (default) with an output vector of length 23 (i.e MOD13Q1 or MYD13Q1) |
nCluster |
Numeric. Number of CPU cores to be used for computational calculations |
outname |
Character vector with the output directory path and filename with extension or only the filename and extension if working directory was set. For example outname = "output_phen.tif". See |
datatype |
Character. Output data type. For example datatype = "INT2S" for signed integer values. See |
rge |
A vector containing minimum and maximum values of the vegetation index used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge = c(0,10000) |
Per pixel, it derives the most recurrent annual Land Surface Phenology (LSP) for a reference period using time series of satellite based greenness values such as the Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). The annual LSP is calculated for all pixels of the input SpatRaster in the same way as for the Phen function. The output is a multiband raster where every band is the expected greenness value at a given time step of the growing season. For example, for MODIS Vegetation Index 16-days composites the number of time steps of the growing season is 23 (frequency = "16-days"), and therefore, the output raster will have 23 bands. A vector with dates for the greenness values is also required
SpatRaster
## Not run: ## DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING## ## Testing with an NDVI spatRast from Central Chile, h = 2## # Load data # SpatRaster library(terra) f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Making the LSP output raster, n bands = 23 # Define the number of cores to be use. In this example we use 1 nc1 <- 1 PhenMap( s = MegaDrought, dates = modis_dates, h = 2, frequency = "16-days", nCluster = nc1, outname = "phen_MD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map1 <- rast("phen_MD.tif")#run only for load phenology brick # plot(map1) ## Testing with an NDVI spatRast from the Atacama Desert, Northern Chile, h=2 ## # Load data SpatRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Dates data("modis_dates") # Making the LSP output raster, n bands = 23 # Define the number of cores to be use. In this example we use 1 nc1 <- 1 PhenMap( s = Bdesert, dates = modis_dates, h = 2, frequency = "16-days", nCluster = 1, outname = "phen_BD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map2 <- rast("phen_BD.tif") #run only for loading the multiband LSP spatRaster # plot(map2) ## End(Not run)## Not run: ## DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING## ## Testing with an NDVI spatRast from Central Chile, h = 2## # Load data # SpatRaster library(terra) f <- system.file("extdata/MegaDrought_spatRast.rda", package = "npphen") MegaDrought <- readRDS(f) # Dates data("modis_dates") # Making the LSP output raster, n bands = 23 # Define the number of cores to be use. In this example we use 1 nc1 <- 1 PhenMap( s = MegaDrought, dates = modis_dates, h = 2, frequency = "16-days", nCluster = nc1, outname = "phen_MD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map1 <- rast("phen_MD.tif")#run only for load phenology brick # plot(map1) ## Testing with an NDVI spatRast from the Atacama Desert, Northern Chile, h=2 ## # Load data SpatRaster f <- system.file("extdata/Bdesert_spatRast.rda", package = "npphen") Bdesert <- readRDS(f) # Dates data("modis_dates") # Making the LSP output raster, n bands = 23 # Define the number of cores to be use. In this example we use 1 nc1 <- 1 PhenMap( s = Bdesert, dates = modis_dates, h = 2, frequency = "16-days", nCluster = 1, outname = "phen_BD.tif", datatype = "INT2S", rge = c(0, 10000) ) # map2 <- rast("phen_BD.tif") #run only for loading the multiband LSP spatRaster # plot(map2) ## End(Not run)
A dataframe with 8-days MODIS NDVI of a pixel with deciduous Nothofagus macrocarpa forest, Central Chile, using MOD13Q1/MYD13Q1 combined
phentsphents
A data frame with 929 rows and 2 variables:
dates for each record
NDVI value for a pixel of deciduous N. macrocarpa forest