PhenAnoMap: PhenAnoMap

Description

Calculates anomalies with respect to the regular phenological cycle using time series (raster) of vegetation greenness.

Usage

PhenAnoMap(s, dates, h, refp, anop, nCluster, outname, format, datatype, rge)

Arguments

Raster stack with greenness (e,g. NDVI or EV) values

dates

Vector with dates at which the greenness values were recorded

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)

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 <U+2013> 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 <U+2013> 19/12/2001). anop y refp can be overlapped

nCluster

Numeric. Number of CPU cores to be used for computational calculations

outname

Character vector with the output path and filename with extension or only the filename and extension if work directory was set. For example outname="output_phen.tif". See writeRaster

format

Character. Output file type. See writeFormats

datatype

Character. Output data type. See dataType

rge

A vector containing minimum and maximum values of the response variable 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)

Value

RasterStack

Details

Similar to PhenAnoma, 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 greeness into 2: the reference period, from which the annual phenological cycle is calculated (same as the Phen 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. It delivers a raster stack with anomalies per date.

Examples

Run this code

# NOT RUN {
##DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING##

## Testing North Hemisphere data. Raster data from Slovenia (EVI index), h=1 ##

# Load data
sl.path<-system.file("extdata/HN_slovenia",package="npphen")
sl_rasters<-list.files(path=sl.path, pattern=glob2rx("slovenia*.tif"), full.names=TRUE)
Slovenia_rasters<-stack(sl_rasters)
sl_dates<-read.csv(system.file("extdata/date_tables/Slovenia_dates.csv", package="npphen"))
Slovenia_dates <- as.Date(sl_dates$date, format='%d/%m/%Y')

# Calculating the anomalies for the last growing season [343:365], refp [1:342]
# In this case refp and anop do not overlap

library(snow)

# Define the number of cores to be use. In this example we use 1
nc1<-1

PhenAnoMap(s=Slovenia_rasters,dates=Slovenia_dates,h=1,refp=c(1:342), anop=c(343:365),
nCluster=nc1,outname="ano_slov.tif", format="GTiff", datatype="FLT4S", rge=c(0,10000))
#map_an1<-raster("ano_slov.tif")#run only for load anomaly brick
#plot(map_an1)


## Testing South Hemisphere data. Raster data from Chile (EVI index), h=2 ##

# Load data
ay.path<-system.file("extdata/HS_aysen",package="npphen")
ayrasters<-list.files(path=ay.path, pattern=glob2rx("aysen*.tif"), full.names=TRUE)
Aysen_rasters<-stack(ayrasters)
ay_dates<-read.csv(system.file("extdata/date_tables/Aysen_dates.csv", package="npphen"))
Aysen_dates <- as.Date(ay_dates$date, format='%d/%m/%Y')

# Making the LSP raster, n bands = 23
# Define the number of cores to be use. In this example we use 1
nc1<-1

PhenAnoMap(s=Aysen_rasters,dates=Aysen_dates,h=2,refp=c(1:354), anop=c(309:331),
nCluster=nc1,outname="ano_aysen.tif", format="GTiff", datatype="FLT4S",rge=c(0,10000))
#map_an2<-raster("ano_aysen.tif")#run only for load anomaly brick
#plot(map_an2)

# }