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
The following code block executes when the notebook is opened.
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/
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")
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
}
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)
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)
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)
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)
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"))
# 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"))
# 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"))
# 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"))
Histograms show the frequency of scences that are available to derive mean values.
To plot additional histograms, change the aoi
to another example area.
Create sets of time series plots for the Julian dates
# NDVI Time Series
These plot use the calendar dates to create NDVI time series ### Okavango Thornscrub
These plot use the calendar dates to create EVI time series ### Okavango Thornscrub
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