compare_raster_time: Compare Temporal Consistency of two Rasters

In antoinestevens/copernicus: Download and Extract COPERNICUS Global Land Vegetation Products

Description

Analysis of the temporal consistency (agreement) of two Raster-class objects with dimensions x (lon), y (lat), z (time). The analysis is performed pixel by pixel. The temporal profile (z dimension) of each pixel at a given (x,y) is extracted in the two raster and compared using a set of statistics. This results in a raster map for each computed statistic, depicting the spatial variability of the temporal coherence (agreement) between the two datasets.

Usage

compare_raster_time(x,y,stats,FUN = NULL,cl = NULL,filename=rasterTmpFile(),...)

Arguments

x	A Raster-class object. Better be a RasterBrick, for faster computing
y	Another Raster-class object to compare with the first one
stats	A character vector with one or more of the following: 'missing', 'cor', 'ax', 'ay', 'bx', 'by', 'ac','acu', 'acs', 'mbe', 'rmsd', 'rmspd', 'rmpdu', 'rmpds','mpdpu','mpdps', 'smoothness','atime','btime'. See details
FUN	An optional function to be applied on the z slot of the input x and y (can be set with setZ, to group statistics according to the resulting classes. If provided, only one stats is allowed. See Examples section.
cl	cluster object for parallel processing. Default is NULL
filename	name of the file where the output raster should be saved. Default is created through rasterTmpFile
...	arguments passed to writeRaster

Details

The stats arguments can take one or more of the following values (see Ji and Gallo, 2006):

• missing Number of missing values in the x and y raster

• cor Correlation coefficient

• ax Intercept of the Geometric Mean Functional Relationship (GMFR) model x = a + by

• ay Intercept of the GMFR model y = a + bx

• bx Coefficient of the GMFR model x = a + by

• by Coefficient of the GMFR model y = a + bx

• ac Agreement Coefficient

• acu Unsystematic Agreement Coefficient

• acs Systematic Agreement Coefficient

• mbe Mean Bias Error

• rmsd Root Mean Square Deviation

• rmspd Root Mean Square Deviation in percentage of mean y

• rmpdu Unsystematic Square Root of Mean Product-difference

• rmpds Systematic Square Root of Mean Product-difference

• mpdpu Proportion of the unsystematic differences over the total difference

• mpdps Proportion of the systematic differences over the total difference

• smoothness Mean absolute value of the first derivative over time

• atime Intercept of the regression over time

• btime Coefficient of the regression over time

Value

A Raster-class object with layers correponding to stats. If FUN argument is provided, layers correponds to the groups resulting from applying FUN to the z slot of input rasters x and y.

Author(s)

Antoine Stevens

References

Ji, L., and Gallo, K. (2006). An Agreement Coefficient for Image Comparison. Photogrammetric Engineering & Remote Sensing 72, 823-833. Meroni, M., Atzberger, C., Vancutsem, C., Gobron, N., Baret, F., Lacaze, R., Eerens, H., and Leo, O. (2013). Evaluation of Agreement Between Space Remote Sensing SPOT-VEGETATION fAPAR Time Series. IEEE Transactions on Geoscience and Remote Sensing 51, 1951-1962. Meroni M., Fasbender D., Balaghi et al. (2015). Testing VGT data continuity between SPOT and PROBA-V missions for operational yield forecasting in North African countries. JRC Technical Report, 28 p.

Examples

## Not run: 
# Let's compare VGT and PROBA-V instruments during their overlaping period
# (October 2013 - March 2014)
# Don't forget to provide in copernicus_options() your user and password details
# for COPERNICUS data portal before running this
# e.g. : copernicus_options(user = "Smith", password = "hello")
# First, get data for NDVI_1km_V1
fn_SPOT <- download_copernicus(product = 'NDVI_1km_V1', begin = '2013-10-01', end = '2014-03-31',
                     tileH = 19, tileV = 4)
# and NDVI_1km_V2 ...
fn_PROBA <- download_copernicus(product = 'NDVI_1km_V2', begin = '2013-10-01', end = '2014-03-31',
                     tileH = 19, tileV = 4)
# Extract NDVI, export to tif
f_SPOT <-extract_copernicus(fn_SPOT,job = "product_comparison",layers = 2)
f_PROBA <-extract_copernicus(fn_PROBA,job = "product_comparison",layers = 1)

# Convert to  rasterBrick
f_SPOT <- sub('\\.h5','_NDVI.tif',f_SPOT)
f_PROBA <- sub('\\.h5','_NDVI.tif',f_PROBA)
SPOT  <- writeRaster(stack(f_SPOT),filename = rasterTmpFile())
PROBA  <- writeRaster(stack(f_PROBA),filename = rasterTmpFile())

# Compare temporal consistency
# Use the 'cl' option to speed up processing!
ct <- compare_raster_time(SPOT,PROBA)
# plot
brks <- c(0,.05,.1,.2,.3,.7)
nb <- length(brks)-1
cols <- rev(heat.colors(nb))
plot(ct,"rmsd",col = cols,breaks = brks) # rmsd
# focus on a small area
# extent corresponding to the Delta of the Po
e <- extent(c(10.5,12.5,44,46)) # xmin,xmax,ymin,ymax
# convert geo coordinates to EPSG 32662 (Platte Carree)
e <- extent(projectExtent(raster(e,crs = crs("+init=epsg:4326")),crs("+init=epsg:32662")))
plot(ct,"ac",ext = e,zlim=c(0,1))

# Let's see if there are different patterns
# according to seasons
# First, define a z slot
d <- scan_file_copernicus(names(SPOT))$Date # get dates
SPOT <- setZ(SPOT,d)
PROBA <- setZ(PROBA,d)
# Compare, grouping by Quarter
# create a special function, appending 'Q' to the quarter number
quarter <- function(x) paste0("Q",lubridate::quarter(x))
quarter(d)
ct_Quarter <- compare_raster_time(SPOT,PROBA,stats = "cor",FUN = quarter)
plot(ct_Quarter)

## End(Not run)

antoinestevens/copernicus documentation built on May 10, 2019, 12:23 p.m.