Setup

Linux Dependencies

If you are running this on a VM outside CyVerse, you may need to install additional dependencies

# sudo add-apt-repository ppa:ubuntugis
# sudo apt-get update
# sudo apt-get install libgdal-dev libgeos-dev libproj-dev
# sudo apt-get install libudunits2-dev

Install missing R Libraries

The following code block executes when the notebook is opened.

Transfer datasets from CyVerse

This command uses wget to download data into the container.

The CyVerse iRODS data store also supports WebDav; we are hosting the extracted raster data in a folder with anonymous public read-only settings.

Data are ~22.8 GB in size, and will take ~18 minutes to transfer within the CyVerse workbench.

# uncomment to run
# cd ~/emsi/data/
# time wget -r -nH --cut-dirs=6 --no-parent -l8 --reject="index.html*" https://data.cyverse.org/dav-anon/iplant/home/tswetnam/emsi/data/collections/

Create Multi-Plot ability

Import Example Data

In the Git repository /data/examples folder I’ve saved data exported from the Google Earth Engine (GEE) NDVI Explorer (.js files provided in the /javascripts folder).

I have extracted the time series of 16-day Landsat 5,7,8 and 5-day Sentinel-2 NDVI data using their Surface Reflectance model with bit-mask flags for cloud, cloud shadow, water, or snow (Landsat), scenes with cloud <20% (Sentinel-2). Even with these filters turned on there are some pixels affected by snow or cloud contaminated pixels.

## Import Landsat 5,7,8 16-day and Sentinel-2 5 day NDVI from Google Earth Engine extracted CSVs

## Okavango Thornscrub
aoi1 <- read.csv("~/emsi/data/examples/l5-okavango-shrubland-ndvi.csv")
names(aoi1) <- c("date", "ndvi")
aoi2 <- read.csv("~/emsi/data/examples/l7-okavango-shrubland-ndvi.csv")
names(aoi2) <- c("date", "ndvi")
aoi3 <- read.csv("~/emsi/data/examples/l8-okavango-shrubland-ndvi.csv")
names(aoi3) <- c("date", "ndvi")
aoi4 <- read.csv("~/emsi/data/examples/s2-okavango-shrubland-ndvi.csv")
names(aoi4) <- c("date", "ndvi")

## Okavango Swamp
aoi5 <- read.csv("~/emsi/data/examples/l5-okavango-swamp-ndvi.csv")
names(aoi5) <- c("date", "ndvi")
aoi6 <- read.csv("~/emsi/data/examples/l7-okavango-swamp-ndvi.csv")
names(aoi6) <- c("date", "ndvi")
aoi7 <- read.csv("~/emsi/data/examples/l8-okavanog-swamp-ndvi.csv")
names(aoi7) <- c("date", "ndvi")
aoi8 <- read.csv("~/emsi/data/examples/s2-okavango-swamp-ndvi.csv")
names(aoi8) <- c("date", "ndvi")

## US Mexico Grassland Mexico
aoi9 <- read.csv("~/emsi/data/examples/l5-usmex-grassmx-ndvi.csv")
names(aoi9) <- c("date", "ndvi")
aoi10 <- read.csv("~/emsi/data/examples/l7-usmex-grassmx-ndvi.csv")
names(aoi10) <- c("date", "ndvi")
aoi11 <- read.csv("~/emsi/data/examples/l8-usmex-grassmx-ndvi.csv")
names(aoi11) <- c("date", "ndvi")
aoi12 <- read.csv("~/emsi/data/examples/s2-usmex-grassmx-ndvi.csv")
names(aoi12) <- c("date", "ndvi")

## US Mexico Grassland USA
aoi13 <- read.csv("~/emsi/data/examples/l5-usmex-grassus-ndvi.csv")
names(aoi13) <- c("date", "ndvi")
aoi14 <- read.csv("~/emsi/data/examples/l7-usmex-grassus-ndvi.csv")
names(aoi14) <- c("date", "ndvi")
aoi15 <- read.csv("~/emsi/data/examples/l8-usmex-grassus-ndvi.csv")
names(aoi15) <- c("date", "ndvi")
aoi16 <- read.csv("~/emsi/data/examples/s2-usmex-grassus-ndvi.csv")
names(aoi16) <- c("date", "ndvi")

## Acre Brazil Agricultural Areas
aoi17 <- read.csv("~/emsi/data/examples/l5-acre-ag-ndvi.csv")
names(aoi17) <- c("date", "ndvi")
aoi18 <- read.csv("~/emsi/data/examples/l7-acre-ag-ndvi.csv")
names(aoi18) <- c("date", "ndvi")
aoi19 <- read.csv("~/emsi/data/examples/l8-acre-ag-ndvi.csv")
names(aoi19) <- c("date", "ndvi")
aoi20 <- read.csv("~/emsi/data/examples/s2-acre-ag-ndvi.csv")
names(aoi20) <- c("date", "ndvi")

## Acre Brazil Rainforest
aoi21 <- read.csv("~/emsi/data/examples/l5-acre-rainforest-ndvi.csv")
names(aoi21) <- c("date", "ndvi")
aoi22 <- read.csv("~/emsi/data/examples/l7-acre-rainforest-ndvi.csv")
names(aoi22) <- c("date", "ndvi")
aoi23 <- read.csv("~/emsi/data/examples/l8-acre-rainforest-ndvi.csv")
names(aoi23) <- c("date", "ndvi")
aoi24 <- read.csv("~/emsi/data/examples/s2-acre-rainforest-ndvi.csv")
names(aoi24) <- c("date", "ndvi")

## Yakutia Taiga
aoi25 <- read.csv("~/emsi/data/examples/l5-yakutia-taiga-ndvi.csv")
names(aoi25) <- c("date", "ndvi")
aoi26 <- read.csv("~/emsi/data/examples/l7-yakutia-taiga-ndvi.csv")
names(aoi26) <- c("date", "ndvi")
aoi27 <- read.csv("~/emsi/data/examples/l8-yakutia-taiga-ndvi.csv")
names(aoi27) <- c("date", "ndvi")
aoi28 <- read.csv("~/emsi/data/examples/s2-yakutia-taiga-ndvi.csv")
names(aoi28) <- c("date", "ndvi")

## Yakutia Alaas
aoi29 <- read.csv("~/emsi/data/examples/l5-yakutia-alaas-ndvi.csv")
names(aoi29) <- c("date", "ndvi")
aoi30 <- read.csv("~/emsi/data/examples/l7-yakutia-alaas-ndvi.csv")
names(aoi30) <- c("date", "ndvi")
aoi31 <- read.csv("~/emsi/data/examples/l8-yakutia-alaas-ndvi.csv")
names(aoi31) <- c("date", "ndvi")
aoi32 <- read.csv("~/emsi/data/examples/s2-yakutia-alaas-ndvi.csv")
names(aoi32) <- c("date", "ndvi")

Calculate Julian dates from calendar

Create a function to set the Julian (1-366) dates from the data, also flags any dates with no data.

# Create Functions for calculating Julian Dates 
aoi.date_julian <-  function(df){
    df$ndvi_nans <- as.numeric(as.character(df$ndvi))
    df$ndvi_range <- ifelse(df$ndvi_nans>0.05,df$ndvi_nans,NA)
    df$ndvi_range <- ifelse(df$ndvi_nans<0.95,df$ndvi_nans,NA)
    df$asdate <- as.Date(df$date, format = "%b %d, %Y")
    df$julian <- yday(df$asdate)
    df$julian_rounded <- round((df$julian/365)*52)*7
    df
}

Execute Functions recursively through all of the aoi

Runs the Julian function on the 32 example data sets.

aoi1 <- aoi.date_julian(aoi1)
aoi2 <- aoi.date_julian(aoi2)
aoi3 <- aoi.date_julian(aoi3)
aoi4 <- aoi.date_julian(aoi4)
aoi5 <- aoi.date_julian(aoi5)
aoi6 <- aoi.date_julian(aoi6)
aoi7 <- aoi.date_julian(aoi7)
aoi8 <- aoi.date_julian(aoi8)
aoi9 <- aoi.date_julian(aoi9)
aoi10 <- aoi.date_julian(aoi10)
aoi11 <- aoi.date_julian(aoi11)
aoi12 <- aoi.date_julian(aoi12)
aoi13 <- aoi.date_julian(aoi13)
aoi14 <- aoi.date_julian(aoi14)
aoi15 <- aoi.date_julian(aoi15)
aoi16 <- aoi.date_julian(aoi16)
aoi17 <- aoi.date_julian(aoi17)
aoi18 <- aoi.date_julian(aoi18)
aoi19 <- aoi.date_julian(aoi19)
aoi20 <- aoi.date_julian(aoi20)
aoi21 <- aoi.date_julian(aoi21)
aoi22 <- aoi.date_julian(aoi22)
aoi23 <- aoi.date_julian(aoi23)
aoi24 <- aoi.date_julian(aoi24)
aoi25 <- aoi.date_julian(aoi25)
aoi26 <- aoi.date_julian(aoi26)
aoi27 <- aoi.date_julian(aoi27)
aoi28 <- aoi.date_julian(aoi28)
aoi29 <- aoi.date_julian(aoi29)
aoi30 <- aoi.date_julian(aoi30)
aoi31 <- aoi.date_julian(aoi31)
aoi32 <- aoi.date_julian(aoi32)

Calculate the loess smoothed mean NDVI for each Julian date and EMSI

Create a function that calculates the EMSI from the NDVI Julian date time series.

# Calculate EMSI for each area of interest
aoi.emsi <-  function(df){
    mean_loess <- predict(loess(ndvi_range ~ julian_rounded, df), df$julian_rounded)
    df$emsi = (df$ndvi_range - mean_loess) / sd.multiperiod(df$ndvi_range,scale=1)
    df
}

Runs the EMSI function on the 32 example data sets.

aoi1 <- aoi.emsi(aoi1)
aoi2 <- aoi.emsi(aoi2)
aoi3 <- aoi.emsi(aoi3)
aoi4 <- aoi.emsi(aoi4)
aoi5 <- aoi.emsi(aoi5)
aoi6 <- aoi.emsi(aoi6)
aoi7 <- aoi.emsi(aoi7)
aoi8 <- aoi.emsi(aoi8)
aoi9 <- aoi.emsi(aoi9)
aoi10 <- aoi.emsi(aoi10)
aoi11 <- aoi.emsi(aoi11)
aoi12 <- aoi.emsi(aoi12)
aoi13 <- aoi.emsi(aoi13)
aoi14 <- aoi.emsi(aoi14)
aoi15 <- aoi.emsi(aoi15)
aoi16 <- aoi.emsi(aoi16)
aoi17 <- aoi.emsi(aoi17)
aoi18 <- aoi.emsi(aoi18)
aoi19 <- aoi.emsi(aoi19)
aoi20 <- aoi.emsi(aoi20)
aoi21 <- aoi.emsi(aoi21)
aoi22 <- aoi.emsi(aoi22)
aoi23 <- aoi.emsi(aoi23)
aoi24 <- aoi.emsi(aoi24)
aoi25 <- aoi.emsi(aoi25)
aoi26 <- aoi.emsi(aoi26)
aoi27 <- aoi.emsi(aoi27)
aoi28 <- aoi.emsi(aoi28)
aoi29 <- aoi.emsi(aoi29)
aoi30 <- aoi.emsi(aoi30)
aoi31 <- aoi.emsi(aoi31)
aoi32 <- aoi.emsi(aoi32)

Merge Landsat 5, 7, 8, and Sentinel-2 time series

This code block creates a single data.frame with all of the NDVI time series. Currently, it is commented out to limit the effect of Sentinel-2 cloud artifacts in the mean NDVI calculation.

#all1 <- rbind(aoi1,aoi2,aoi3,aoi4)
#all2 <- rbind(aoi5,aoi6,aoi7,aoi8)
#all3 <- rbind(aoi9,aoi10,aoi11,aoi12)
#all4 <- rbind(aoi13,aoi14,aoi15,aoi16)
#all5 <- rbind(aoi17,aoi18,aoi19,aoi20)
#all6 <- rbind(aoi21,aoi22,aoi23,aoi24)
#all7 <- rbind(aoi25,aoi26,aoi27,aoi28)
#all8 <- rbind(aoi29,aoi30,aoi31,aoi32)

Merge Landsat 5, 7, 8 time series

all1 <- rbind(aoi1,aoi2,aoi3)
all2 <- rbind(aoi5,aoi6,aoi7)
all3 <- rbind(aoi9,aoi10,aoi11)
all4 <- rbind(aoi13,aoi14,aoi15)
all5 <- rbind(aoi17,aoi18,aoi19)
all6 <- rbind(aoi21,aoi22,aoi23)
all7 <- rbind(aoi25,aoi26,aoi27)
all8 <- rbind(aoi29,aoi30,aoi31)

Calculate the Loess regression spine for the mean NDVI and

We use a loess regression to show the running mean for the NDVI and EMSI time series. ### Okavango

# Calculate mean NDVI for Julian period in Shrubland
ndvi_mean_shrubland <- c(NA)
ndvi_mean_shrubland$l5 <- setNames(aggregate(aoi1$ndvi_range, list(aoi1$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_shrubland$l7 <- setNames(aggregate(aoi2$ndvi_range, list(aoi2$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_shrubland$l8 <- setNames(aggregate(aoi3$ndvi_range, list(aoi3$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_shrubland$s2 <- setNames(aggregate(aoi4$ndvi_range, list(aoi4$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_shrubland$all <- setNames(aggregate(all1$ndvi_range, list(all1$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
# Calculate mean NDVI for Julian period in Swamp
ndvi_mean_swamp <- c(NA)
ndvi_mean_swamp$l5 <- setNames(aggregate(aoi5$ndvi_range, list(aoi5$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_swamp$l7 <- setNames(aggregate(aoi6$ndvi_range, list(aoi6$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_swamp$l8 <- setNames(aggregate(aoi7$ndvi_range, list(aoi7$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_swamp$s2 <- setNames(aggregate(aoi8$ndvi_range, list(aoi8$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_swamp$all <- setNames(aggregate(all2$ndvi_range, list(all2$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))

US Mexico

# Calculate mean NDVI for Julian period in Mexico
ndvi_mean_grassmx <- c(NA)
ndvi_mean_grassmx$l5 <- setNames(aggregate(aoi9$ndvi_range, list(aoi9$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassmx$l7 <- setNames(aggregate(aoi10$ndvi_range, list(aoi10$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassmx$l8 <- setNames(aggregate(aoi11$ndvi_range, list(aoi11$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassmx$s2 <- setNames(aggregate(aoi12$ndvi_range, list(aoi12$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassmx$all <- setNames(aggregate(all3$ndvi_range, list(all3$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
# Calculate mean NDVI for Julian period in USA
ndvi_mean_grassus <- c(NA)
ndvi_mean_grassus$l5 <- setNames(aggregate(aoi13$ndvi_range, list(aoi13$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassus$l7 <- setNames(aggregate(aoi14$ndvi_range, list(aoi14$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassus$l8 <- setNames(aggregate(aoi15$ndvi_range, list(aoi15$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassus$s2 <- setNames(aggregate(aoi16$ndvi_range, list(aoi16$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_grassus$all <- setNames(aggregate(all4$ndvi_range, list(all4$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))

Acre

# Calculate mean NDVI for Julian period in Agricultural
ndvi_mean_ag <- c(NA)
ndvi_mean_ag$l5 <- setNames(aggregate(aoi17$ndvi_range, list(aoi17$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_ag$l7 <- setNames(aggregate(aoi18$ndvi_range, list(aoi18$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_ag$l8 <- setNames(aggregate(aoi19$ndvi_range, list(aoi19$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_ag$s2 <- setNames(aggregate(aoi20$ndvi_range, list(aoi20$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_ag$all <- setNames(aggregate(all5$ndvi_range, list(all5$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
# Calculate mean NDVI for Julian period in Rainforest
ndvi_mean_rainforest <- c(NA)
ndvi_mean_rainforest$l5 <- setNames(aggregate(aoi21$ndvi_range, list(aoi21$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_rainforest$l7 <- setNames(aggregate(aoi22$ndvi_range, list(aoi22$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_rainforest$l8 <- setNames(aggregate(aoi23$ndvi_range, list(aoi23$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_rainforest$s2 <- setNames(aggregate(aoi24$ndvi_range, list(aoi24$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_rainforest$all <- setNames(aggregate(all6$ndvi_range, list(all6$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))

Yakutia

# Calculate mean NDVI for Julian period in Taiga
ndvi_mean_taiga <- c(NA)
ndvi_mean_taiga$l5 <- setNames(aggregate(aoi25$ndvi_range, list(aoi25$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_taiga$l7 <- setNames(aggregate(aoi26$ndvi_range, list(aoi26$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_taiga$l8 <- setNames(aggregate(aoi27$ndvi_range, list(aoi27$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_taiga$s2 <- setNames(aggregate(aoi28$ndvi_range, list(aoi28$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_taiga$all <- setNames(aggregate(all7$ndvi_range, list(all7$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
# Calculate mean NDVI for Julian period in Alaas
ndvi_mean_alaas <- c(NA)
ndvi_mean_alaas$l5 <- setNames(aggregate(aoi29$ndvi_range, list(aoi29$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_alaas$l7 <- setNames(aggregate(aoi30$ndvi_range, list(aoi30$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_alaas$l8 <- setNames(aggregate(aoi31$ndvi_range, list(aoi31$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_alaas$s2 <- setNames(aggregate(aoi32$ndvi_range, list(aoi32$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))
ndvi_mean_alaas$all <- setNames(aggregate(all8$ndvi_range, list(all8$julian_rounded), mean, na.rm=TRUE, na.action=NULL), c("julian", "ndvi_mean"))

Plot Scene Frequency Histograms

Histograms show the frequency of scences that are available to derive mean values.

To plot additional histograms, change the aoi to another example area.

Plot Time Series

Create sets of time series plots for the Julian dates

Okavango Thornscrub with Legend

### Okavango Thornscrub

Okavango Swamp

US-Mex Border Mexico

### US-Mex Border USA
### Acre Agriculture
### Acre Rainforest
### Yakutia Taiga
### Yakutia Alaas
### Okavango Thornscrub

Okavango Swamp

Mexico Grassland

USA Grassland

Acre Agriculture

Acre Rainforest

Yakutia Taiga

Yakutia Alaas

Plot all Julian NDVI and EMSI

# NDVI Time Series

These plot use the calendar dates to create NDVI time series ### Okavango Thornscrub

### Okavango Swamp
### Mexico Grassland
### USA Grassland
### Acre Agriculture
### Acre Rainforest
### Yakutia Taiga
### Yakutia Alaas

EMSI Time Series

These plot use the calendar dates to create EVI time series ### Okavango Thornscrub

### Okavango Swamp
### Mexico Grassland
### USA Grassland
### Acre Agriculture
### Acre Rainforest
### Yakutia Taiga
### Yakutia Alaas
### Plot all Calendar NDVI and EMSI

m1 <- multiplot(emsi_shrubland, ndvi_shrubland, emsi_swamp, ndvi_swamp, emsi_grassmx, ndvi_grassmx, emsi_grassus, ndvi_grassus, cols=2)
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 418
## Warning: Removed 5 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 415
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1 rows containing missing values (geom_point).

m2 <- multiplot(emsi_acreag, ndvi_acreag, emsi_acrerainforest, ndvi_acrerainforest, emsi_taiga, ndvi_taiga, emsi_alaas, ndvi_alaas, cols=2)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 21 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 463
## Warning: Removed 13 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 476
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 17327
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 8
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.0191e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 312
## Warning: Removed 1 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 17327
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 8
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.0191e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 313
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 300
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 300

m1 <- multiplot(emsi_shrubland, ndvi_shrubland,cols=1)
## `geom_smooth()` using formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : k-d tree limited by memory. ncmax= 418