4. An environmental index with a Principal Component Analysis

Jan 3, 2023 · 3 min read · pca ·

The objective of this demo is to create a new environment index from the analysis of the principal components of multiple variables in the city of San Diego, the dataset used is:

Landcover
% Impervious
% Tree cover
NDVI index
EVI index

All the dataset is from USGS National Land Cover Database and MOD13Q1.006.

Information:

- For this demo you need to have rgee, sf, tidyverse, and viridis packages previously installed.

4.1 Requirements

library(rgee)
library(tidyverse)
library(sf)
library(viridis)
source("https://raw.githubusercontent.com/ambarja/rgee_landcover/main/pca_rgee.R")
ee_Initialize()

── rgee 1.1.2 ────────────────────────── earthengine-api 0.1.292 ── 
 ✓ user: not_defined
 ✓ Initializing Google Earth Engine:  DONE!
 ✓ Earth Engine account: users/ambarja 
───────────────────────────────────────────────────────────────────

4.2 Vector layer reading of our study area

geodata <- st_read(
  "https://github.com/healthinnovation/sdb-gpkg/raw/main/SanDiego.gpkg"
  )
sandiego <- geodata %>% sf_as_ee()

4.3 Processing data with rgee

years <- c(1992, 2001,2004,2006,2008,2011,2013,2016) %>% ee$List()
nldc <- years$
  map(ee_utils_pyfunc(function(x) {
    ee$ImageCollection("USGS/NLCD_RELEASES/2016_REL")$
      select(c("landcover"))$
      filter(ee$Filter$calendarRange(x, x, "year"))$
      sum()$
      clip(sandiego)
  })
  )

nldc <- ee$ImageCollection(nldc)$
  toBands()$
  select("7_landcover")$
  rename("nldc")

years <- c(2011, 2016) %>% ee$List()
impervious <- years$
  map(ee_utils_pyfunc(function(x) {
    ee$ImageCollection("USGS/NLCD_RELEASES/2016_REL")$
      select("impervious")$
      filter(ee$Filter$calendarRange(x, x, "year"))$
      sum()$
      clip(sandiego)
  })
  )

imp <- ee$ImageCollection(impervious)$
  toBands()$
  select("0_impervious")$
  rename("imp")

years <- c(2011, 2016) %>% ee$List()

tree_cover <- years$
  map(ee_utils_pyfunc(function(x) {
    ee$ImageCollection("USGS/NLCD_RELEASES/2016_REL")$
      select("percent_tree_cover")$
      filter(ee$Filter$calendarRange(x, x, "year"))$
      sum()$
      clip(sandiego)
  })
  )

tree <- ee$ImageCollection(tree_cover)$
  toBands()$
  select("1_percent_tree_cover")$
  rename("tree")

modis_years <- c(2010:2021) %>% ee$List()
ndvi <- modis_years$
  map(ee_utils_pyfunc(function(x){
    ee$ImageCollection("MODIS/006/MOD13Q1")$
      select("NDVI")$
      filter(ee$Filter$calendarRange(x,x,"year"))$
      mean()$
      multiply(0.0001)$
      clip(sandiego)
  })
  )

ndvi <- ee$ImageCollection(ndvi)$
  toBands()$
  select("0_NDVI")$
  rename("ndvi") 

evi <-  modis_years$
  map(ee_utils_pyfunc(function(x){
    ee$ImageCollection("MODIS/006/MOD13Q1")$
      select("EVI")$
      filter(ee$Filter$calendarRange(x,x,"year"))$
      mean()$
      multiply(0.0001)$
      clip(sandiego)
  })
  )

evi <- ee$ImageCollection(evi)$
  toBands()$
  select("0_EVI")$
  rename("evi") 

stack_evironment <- ee$Image(
    c(nldc, imp, tree, ndvi, evi)
    )$
    toDouble()

4.4 Variable standardization

standCov <-
  ee_scale(
    image = stack_evironment,
    ee_feature = sandiego,
    scale = 30,
    namevar = c("nldc","imp","tree","ndvi","evi")
  )

4.5 PCA

PCA <-
  ee_pca(
    image = standCov,
    ee_feature = sandiego,
    scale = 30,
    nvar = 5
  )$select(sprintf("pc%1$s", 1:5))
Map$centerObject(sandiego)
Map$addLayer(PCA,visParams = list(bands=c("pc1","pc2","pc3")))

4.6 Table of eingvectors

eVectors(
  image = standCov,
  ee_feature = sandiego,
  scale = 30,
  nvar = 5
)

  eingvector  eVec1  eVec2  eVec3  eVec4  eVec5
1          1  0.342 -0.237  0.354  0.597  0.588
2          2  0.598 -0.661 -0.355 -0.193 -0.205
3          3 -0.141 -0.395  0.815 -0.257 -0.306
4          4  0.711  0.592  0.290 -0.177 -0.170
5          5  0.001  0.024 -0.026  0.714 -0.700

4.7 Table of eingvalues

eValues(
  image = standCov,
  ee_feature = sandiego,
  scale = 30,
  nvar = 5
)

  eingvalue values
1     eVal1   2.40
2     eVal2   1.38
3     eVal3   0.74
4     eVal4   0.42
5     eVal5   0.06

4.8 Summary table

table <-
  imporPCA(
    image = standCov,
    ee_feature = sandiego,
    scale = 30,
    nvar = 5
  )
table

  eingvalue values component variance cumulative
1     eVal1   2.40       pc1     48.0       48.0
2     eVal2   1.38       pc2     27.6       75.6
3     eVal3   0.74       pc3     14.8       90.4
4     eVal4   0.42       pc4      8.4       98.8
5     eVal5   0.06       pc5      1.2      100.0

4.9 Plot PCA

table %>%
  ggplot(aes(x = reorder(component, variance), y = variance, fill = variance)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  coord_flip() +
  labs(x = "PCs", y = "%Total variance")

4.10 Construction of an environmental index

index <- PCA$
  select(
  c("pc1","pc2","pc3")
  )$
  reduce("sum")

4.11 Identifying min and max values of environmental index

(minmax <- index$reduceRegion(
  reducer = ee$Reducer$minMax(),
  geometry = sandiego,
  scale = 50*1000)$
  getInfo())

$sum_max
[1] 0.9375844

$sum_min
[1] -0.6983301

4.12 Mapping environmental index

viz <- list(
  min = -0.6983301 ,
  max = 0.9375844,
  palette = viridis(n = 100)
)

Map$addLayer(index,visParams = viz) + 
  Map$addLegend(visParams = viz,name = "Environment index")