Alaska
NDVI analysis {remotePARTS}
Clay Morrow
-2022-06-27
+2023-01-20
+Introduction
needed:devtools::install_github("morrowcj/remotePARTS")Then, ensure that the package is loaded into your library:
-library(remotePARTS)library(remotePARTS)This vignette will use dplyr and ggplot2
for visualizing the data:
library(dplyr)
-##
-## Attaching package: 'dplyr'
-## The following objects are masked from 'package:stats':
-##
-## filter, lag
-## The following objects are masked from 'package:base':
-##
-## intersect, setdiff, setequal, union
-library(ggplot2)
-## Warning: package 'ggplot2' was built under R version 4.0.5library(dplyr)
+##
+## Attaching package: 'dplyr'
+## The following objects are masked from 'package:stats':
+##
+## filter, lag
+## The following objects are masked from 'package:base':
+##
+## intersect, setdiff, setequal, union
+library(ggplot2)Alaska datasets
@@ -374,8 +393,8 @@Alaska datasets
maps).For this vignette, we well also create a smaller 3000 pixel subsample
ndvi_AK3000 for demonstrative purposes:
data("ndvi_AK10000")
-ndvi_AK3000 <- ndvi_AK10000[seq_len(3000),] # first 3000 pixels from the random 10Kdata("ndvi_AK10000")
+ndvi_AK3000 <- ndvi_AK10000[seq_len(3000),] # first 3000 pixels from the random 10Kndvi_AK10000 is a data.frame with 37
columns. lng and lat are longitude and
latitude, respectively. AR_coef and CLS_coef
@@ -391,67 +410,67 @@
Alaska datasets
not handle any missing data. It is essential that all missing data are removed prior to conducting any analyses. -str(ndvi_AK10000)
-## 'data.frame': 10000 obs. of 37 variables:
-## $ lng : num -151 -145 -144 -163 -145 ...
-## $ lat : num 63.3 64.7 62.7 67.7 65.8 ...
-## $ AR_coef : num 0.01071 0.003586 0.002488 0.000373 0.002505 ...
-## $ CLS_coef: num 0.00529 0.001537 0.002276 0.000476 0.001514 ...
-## $ land : Factor w/ 10 levels "Evergr needle",..: 6 7 7 6 7 7 6 6 7 7 ...
-## $ ndvi1982: num 1.88 7.68 8.46 6.85 10.25 ...
-## $ ndvi1983: num 1.74 8.14 7.99 6.97 10.35 ...
-## $ ndvi1984: num 1.96 8.47 8.04 6.67 10.29 ...
-## $ ndvi1985: num 1.74 7.72 7.96 6.9 10.02 ...
-## $ ndvi1986: num 1.92 8.31 8.64 6.43 10.39 ...
-## $ ndvi1987: num 2.04 8.89 8.4 7 10.47 ...
-## $ ndvi1988: num 2.13 8.5 8.32 7.12 10.05 ...
-## $ ndvi1989: num 1.93 8.99 8.79 6.24 10.16 ...
-## $ ndvi1990: num 1.73 8.66 8.45 7.46 10.31 ...
-## $ ndvi1991: num 1.99 8.74 8.34 6.86 10.28 ...
-## $ ndvi1992: num 1.74 7.99 7.94 6.53 9.77 ...
-## $ ndvi1993: num 1.98 8.62 7.93 6.97 10.35 ...
-## $ ndvi1994: num 2.04 8.59 8.26 7.15 10.25 ...
-## $ ndvi1995: num 1.9 8.1 8.55 7.08 10.18 ...
-## $ ndvi1996: num 2.09 8.61 8.23 7.08 10.33 ...
-## $ ndvi1997: num 1.81 8.77 8.99 7.24 10.23 ...
-## $ ndvi1998: num 1.95 7.9 8.57 7.45 9.99 ...
-## $ ndvi1999: num 1.75 8.14 8.61 6.75 10.04 ...
-## $ ndvi2000: num 1.71 7.5 8.37 6.7 9.92 ...
-## $ ndvi2001: num 1.57 7.76 8.29 6.66 9.69 ...
-## $ ndvi2002: num 1.72 7.83 8.62 7.09 9.76 ...
-## $ ndvi2003: num 1.95 7.94 8.35 6.87 9.76 ...
-## $ ndvi2004: num 2.2 8.37 8.83 7.3 11.74 ...
-## $ ndvi2005: num 2.18 8.45 8.63 7.18 11.62 ...
-## $ ndvi2006: num 2.06 8.32 8.36 6.74 11.59 ...
-## $ ndvi2007: num 2.39 8.53 8.63 7.47 11.91 ...
-## $ ndvi2008: num 2.36 8.6 8.41 7 11.8 ...
-## $ ndvi2009: num 2.35 10.02 9.27 6.59 11.62 ...
-## $ ndvi2010: num 2.95 11.09 9.26 7.09 11.46 ...
-## $ ndvi2011: num 2.5 8.9 8.94 7.05 12.15 ...
-## $ ndvi2012: num 2.71 8.71 8.39 6.21 10.97 ...
-## $ ndvi2013: num 2.54 9 8.95 6.93 10.38 ...str(ndvi_AK10000)
+## 'data.frame': 10000 obs. of 37 variables:
+## $ lng : num -151 -145 -144 -163 -145 ...
+## $ lat : num 63.3 64.7 62.7 67.7 65.8 ...
+## $ AR_coef : num 0.01071 0.003586 0.002488 0.000373 0.002505 ...
+## $ CLS_coef: num 0.00529 0.001537 0.002276 0.000476 0.001514 ...
+## $ land : Factor w/ 10 levels "Evergr needle",..: 6 7 7 6 7 7 6 6 7 7 ...
+## $ ndvi1982: num 1.88 7.68 8.46 6.85 10.25 ...
+## $ ndvi1983: num 1.74 8.14 7.99 6.97 10.35 ...
+## $ ndvi1984: num 1.96 8.47 8.04 6.67 10.29 ...
+## $ ndvi1985: num 1.74 7.72 7.96 6.9 10.02 ...
+## $ ndvi1986: num 1.92 8.31 8.64 6.43 10.39 ...
+## $ ndvi1987: num 2.04 8.89 8.4 7 10.47 ...
+## $ ndvi1988: num 2.13 8.5 8.32 7.12 10.05 ...
+## $ ndvi1989: num 1.93 8.99 8.79 6.24 10.16 ...
+## $ ndvi1990: num 1.73 8.66 8.45 7.46 10.31 ...
+## $ ndvi1991: num 1.99 8.74 8.34 6.86 10.28 ...
+## $ ndvi1992: num 1.74 7.99 7.94 6.53 9.77 ...
+## $ ndvi1993: num 1.98 8.62 7.93 6.97 10.35 ...
+## $ ndvi1994: num 2.04 8.59 8.26 7.15 10.25 ...
+## $ ndvi1995: num 1.9 8.1 8.55 7.08 10.18 ...
+## $ ndvi1996: num 2.09 8.61 8.23 7.08 10.33 ...
+## $ ndvi1997: num 1.81 8.77 8.99 7.24 10.23 ...
+## $ ndvi1998: num 1.95 7.9 8.57 7.45 9.99 ...
+## $ ndvi1999: num 1.75 8.14 8.61 6.75 10.04 ...
+## $ ndvi2000: num 1.71 7.5 8.37 6.7 9.92 ...
+## $ ndvi2001: num 1.57 7.76 8.29 6.66 9.69 ...
+## $ ndvi2002: num 1.72 7.83 8.62 7.09 9.76 ...
+## $ ndvi2003: num 1.95 7.94 8.35 6.87 9.76 ...
+## $ ndvi2004: num 2.2 8.37 8.83 7.3 11.74 ...
+## $ ndvi2005: num 2.18 8.45 8.63 7.18 11.62 ...
+## $ ndvi2006: num 2.06 8.32 8.36 6.74 11.59 ...
+## $ ndvi2007: num 2.39 8.53 8.63 7.47 11.91 ...
+## $ ndvi2008: num 2.36 8.6 8.41 7 11.8 ...
+## $ ndvi2009: num 2.35 10.02 9.27 6.59 11.62 ...
+## $ ndvi2010: num 2.95 11.09 9.26 7.09 11.46 ...
+## $ ndvi2011: num 2.5 8.9 8.94 7.05 12.15 ...
+## $ ndvi2012: num 2.71 8.71 8.39 6.21 10.97 ...
+## $ ndvi2013: num 2.54 9 8.95 6.93 10.38 ...For this demonstration, we are interested in asking the following questions using these data: “Is NDVI in Alaska increasing over time?”; “Are Alaska’s NDVI time trends associated with land-cover classes?”; and “Do Alaska’s NDVI time trends differ with latitude?”
The figure below shows a temporal cross-section of these data for 1982, 1998, and 2013.
-reshape2::melt(ndvi_AK10000, measure = c("ndvi1982", "ndvi1998", "ndvi2013")) %>%
- ggplot(aes(x = lng, y = lat, col = value )) +
- geom_point(size = .1) +
- labs(col = "ndvi") +
- facet_wrap(~ gsub("ndvi", "", variable), ncol = 3) +
- scale_color_viridis_c(option = "magma") +
- labs(x = "Longitude", y = "Latitude")
reshape2::melt(ndvi_AK10000, measure = c("ndvi1982", "ndvi1998", "ndvi2013")) %>%
+ ggplot(aes(x = lng, y = lat, col = value )) +
+ geom_point(size = .1) +
+ labs(col = "ndvi") +
+ facet_wrap(~ gsub("ndvi", "", variable), ncol = 3) +
+ scale_color_viridis_c(option = "magma") +
+ labs(x = "Longitude", y = "Latitude")
The following figure shows how Alaska’s three primary land-cover classes are distributed.
-ndvi_AK10000 %>%
-ggplot(aes(x = lng, y = lat, col = land)) +
- geom_point(size = .1) +
- scale_color_viridis_d(direction = -1, end = .9) +
- labs(y = "Latitude", x = "Longitude", col = "Land cover", fill = "Land cover")
ndvi_AK10000 %>%
+ggplot(aes(x = lng, y = lat, col = land)) +
+ geom_point(size = .1) +
+ scale_color_viridis_d(direction = -1, end = .9) +
+ labs(y = "Latitude", x = "Longitude", col = "Land cover", fill = "Land cover")
Use help("ndvi_AK) to see documentation for these
datasets.
Time-series analysis
maximum likelihood (REML). To do so, we must extract only our NDVI columns as the matrixY. We’ll do this by matching all
column names containing “ndvi” and slicing the data.frame:
-ndvi.cols <- grep("ndvi", names(ndvi_AK3000), value = TRUE)
-Y <- as.matrix(ndvi_AK3000[, ndvi.cols])ndvi.cols <- grep("ndvi", names(ndvi_AK3000), value = TRUE)
+Y <- as.matrix(ndvi_AK3000[, ndvi.cols])We also need a 2-column coordinate matrix coords:
coords <- as.matrix(ndvi_AK3000[, c("lng", "lat")])coords <- as.matrix(ndvi_AK3000[, c("lng", "lat")])Y and coords are then passed to
fitAR_map() with default settings:
ARfit <- fitAR_map(Y = Y, coords = coords)ARfit <- fitAR_map(Y = Y, coords = coords)Coefficient estimates can be obtained from ARfit with
coefficients(). The first column is the estimate of \(\beta_0\), \(\hat{\beta_0}\), and the second is \(\hat{\beta_1}\).
head(coefficients(ARfit))
-## (Intercept) t
-## [1,] 1.703446 0.021923615
-## [2,] 7.980618 0.030464243
-## [3,] 8.146784 0.021130056
-## [4,] 6.883592 0.002582853
-## [5,] 10.081968 0.026466877
-## [6,] 7.628119 -0.010395964head(coefficients(ARfit))
+## (Intercept) t
+## [1,] 1.703446 0.021923615
+## [2,] 7.980618 0.030464243
+## [3,] 8.146784 0.021130056
+## [4,] 6.883592 0.002582853
+## [5,] 10.081968 0.026466877
+## [6,] 7.628119 -0.010395964These time-series analyses calculate the time trend in the raw NDVI data. In most situations it makes sense to ask if there are time trends in the relative NDVI values, that is, changes in NDVI relative to the @@ -520,22 +539,22 @@
Time-series analysis
ARfit$coefficients, and since the coefficients for the
trend are in the column of the coefficients matrix named t,
the scaling is performed as
-ARfit$coefficients[, "t"] <- ARfit$coefficients[,"t"]/rowMeans(ndvi_AK3000[, ndvi.cols])
-ndvi_AK3000$AR_coef <- coefficients(ARfit)[, "t"] # save time trend coefficientARfit$coefficients[, "t"] <- ARfit$coefficients[,"t"]/rowMeans(ndvi_AK3000[, ndvi.cols])
+ndvi_AK3000$AR_coef <- coefficients(ARfit)[, "t"] # save time trend coefficientThese scaled values of the time trend are then stored in the
ndvi_AK3000 data frame.
Below is an image of the estimated coefficients (pre-calculated and
scaled) for the full ndvi_AK10000. From this, it appears
that northern latitudes may be greening faster than more southern
latitudes.
ndvi_AK10000 %>%
- ggplot(aes(x = lng, y = lat, col = AR_coef)) +
- geom_point(size = .1) +
- scale_color_gradient2(high = "red", low = "blue",
- mid = "grey90", midpoint = 0) +
- guides(fill = "none") +
- labs(y = "Latitude", x = "Longitude", col = expression(beta[1]))
ndvi_AK10000 %>%
+ ggplot(aes(x = lng, y = lat, col = AR_coef)) +
+ geom_point(size = .1) +
+ scale_color_gradient2(high = "red", low = "blue",
+ mid = "grey90", midpoint = 0) +
+ guides(fill = "none") +
+ labs(y = "Latitude", x = "Longitude", col = expression(beta[1]))
fitAR_map and its conditional least-squares counterpart,
fitCLS_map, are wrappers for the functions
fitAR and fitCLS which conduct individual time
@@ -556,7 +575,7 @@
Distance
The first step is to calculate the distances among pixels as a
distance matrix D. Here, we’ll calculate relative distances
with distm_scaled() from our coordinate matrix.
D <- distm_scaled(coords)D <- distm_scaled(coords)distm_scaled() scales distances across the spatial
domain so that the greatest distance between two pixels is 1. Note that
because distm_scaled() uses
@@ -580,13 +599,13 @@
Covariance
that the covariance structure among pixels is given by \(\Sigma = \eta I + (1-\eta)V\) where \(I\) is the identity matrix.If we know the range parameter \(r\), we can calculate V from
D with covar_exp():
r <- 0.1
-V <- covar_exp(D, r)r <- 0.1
+V <- covar_exp(D, r)And we could add a known nugget to V to obtain
Sigma:
nugget <- 0.2
-I <- diag(nrow(V)) # identity matrix
-Sigma <- nugget*I + (1-nugget)*Vnugget <- 0.2
+I <- diag(nrow(V)) # identity matrix
+Sigma <- nugget*I + (1-nugget)*VSee ?covar_exp() for a description of the covariance
functions provided by remotePARTS and for more information
regarding their use.
Known parameters
V. This means that even with only the 3000-pixel subset, it
will take a few minutes to finish the computations on most
computers.
-GLS.0 <- fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V, nugget = nugget)GLS.0 <- fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V, nugget = nugget)Note also that fitGLS adds the nugget to V
internally. If we wanted to do this ourselves, we could pass the
covariance matrix Sigma which already contains a nugget
component and then set the nugget argument of
fitGLS to 0:
fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = Sigma, nugget = 0) # equivalentfitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = Sigma, nugget = 0) # equivalentThe estimates for our land class effects can be extracted with
coefficients().
coefficients(GLS.0)
-## landShrubland landSavanna landGrassland
-## 0.0008646565 0.0002554862 0.0011294576coefficients(GLS.0)
+## landShrubland landSavanna landGrassland
+## 0.0008646565 0.0002554862 0.0011294576The full model fit is given by
-GLS.0
-##
-## Call:
-## fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V,
-## nugget = nugget)
-##
-## t-tests:
-## Est t.stat pval.t
-## landShrubland 0.0008647 0.5065 0.6126
-## landSavanna 0.0002555 0.1493 0.8813
-## landGrassland 0.0011295 0.6648 0.5062
-##
-## F-test:
-## model df_F SSE MSE logLik Fstat pval_F
-## 1 AR_coef ~ 0 + land 2 0.1128 3.764e-05 12808 4.461 0.01163
-## 2 AR_coef ~ 1 2997 0.1131 3.775e-05 12802 NA NAGLS.0
+##
+## Call:
+## fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V,
+## nugget = nugget)
+##
+## t-tests:
+## Est t.stat pval.t
+## landShrubland 0.0008647 0.5065 0.6126
+## landSavanna 0.0002555 0.1493 0.8813
+## landGrassland 0.0011295 0.6648 0.5062
+##
+## F-test:
+## model df_F SSE MSE logLik Fstat pval_F
+## 1 AR_coef ~ 0 + land 2 0.1128 3.764e-05 12808 4.461 0.01163
+## 2 AR_coef ~ 1 2997 0.1131 3.775e-05 12802 NA NANote that, although none of the three land-cover classes shows a statistically significant time trend in (proportional) NDVI, the F-test shows that there is a statistically significant difference among @@ -687,12 +706,11 @@
Spatial parameters
start the search for the optimal range parameter at 0.1. For this example, we will also specifyfit.n = 3000, which ensures
that all pixels are used to estimate spatial parameters.
-corfit <- fitCor(resids = residuals(ARfit), coords = coords,
- covar_FUN = "covar_exp", start = list(range = 0.1),
- fit.n = 3000)
-(range.opt = corfit$spcor)
-## range
-## 0.1083729corfit <- fitCor(resids = residuals(ARfit), coords = coords, covar_FUN = "covar_exp",
+ start = list(range = 0.1), fit.n = 3000)
+(range.opt = corfit$spcor)
+## range
+## 0.1083729By default, fitCor() uses distm_scaled() to
calculate distances from the coordinate matrix, but any function that
returns a distance matrix can be specified with the
@@ -701,10 +719,9 @@
Spatial parameters
instead usedistm_km() to calculate distance in km instead
of relative distances, we would need to scale our starting range
parameter by the maximum distance in km of our map:
-max.dist <- max(distm_km(coords))
-corfit.km <- fitCor(resids = residuals(ARfit), coords = coords,
- covar_FUN = "covar_exp", start = list(range = max.dist*0.1),
- distm_FUN = "distm_km", fit.n = 3000)max.dist <- max(distm_km(coords))
+corfit.km <- fitCor(resids = residuals(ARfit), coords = coords, covar_FUN = "covar_exp",
+ start = list(range = max.dist*0.1), distm_FUN = "distm_km", fit.n = 3000)Note that, depending on the covariance function used, not all
parameters will need scaling. For example, covar_exppow()
is an exponential-power covariance function and takes a range and shape
@@ -712,7 +729,7 @@
Spatial parameters
?fitCor() for more details.
After we’ve obtained our range parameter estimate, we can use it to
re-calculate the V matrix:
V.opt <- covar_exp(D, range.opt)V.opt <- covar_exp(D, range.opt)