From 4b7c6981e5e6224d034ab3756ac0c2e5acdc21d0 Mon Sep 17 00:00:00 2001 From: Dainius Date: Sat, 19 Mar 2022 15:03:10 +0100 Subject: [PATCH] Update text and start porting to terra --- index.Rmd | 30 +++-- index.html | 375 ++++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 305 insertions(+), 100 deletions(-) diff --git a/index.Rmd b/index.Rmd index c33cba3..a06db19 100644 --- a/index.Rmd +++ b/index.Rmd @@ -30,7 +30,7 @@ a, a:visited { # MODIS based time series analysis using BFAST Monitor and BFAST Lite -![WUR logo](https://www.wur.nl/upload/5e55a1d9-47c3-4149-9384-26d5147ee2d7_WUR_RGB_standard.png) +WUR logo [Jan Verbesselt, Dainius Masiliūnas](http://www.wur.eu/changemonitor) @@ -102,7 +102,7 @@ pkgTest <- function(x) } library(x, character.only = TRUE) } -neededPackages <- c("zoo", "bfast", "raster", "leaflet", "MODISTools") +neededPackages <- c("zoo", "bfast", "terra", "raster", "leaflet", "MODISTools") for (package in neededPackages){pkgTest(package)} ``` @@ -111,9 +111,9 @@ Loading extra function `timeser()` to create a time series object in R: ```{r, message = FALSE} -# Function to create time series object +# Utility function to create time series object from a numeric vector # val_array: data array for one single pixel (length is number of time steps) -# time_array: array with dates at which raster data is recorded (same length as val_array) +# time_array: array with Dates at which raster data is recorded (same length as val_array) timeser <- function(val_array, time_array) { z <- zoo(val_array, time_array) # create zoo object yr <- as.numeric(format(time(z), "%Y")) # extract the year numbers @@ -176,8 +176,8 @@ Second, we create a raster brick using the `mt_to_raster` function that is inclu ```{r mttoraster} # convert df to raster -VI_r <- mt_to_raster(df = VI) -QA_r <- mt_to_raster(df = QA) +VI_r <- rast(mt_to_raster(df = VI)) +QA_r <- rast(mt_to_raster(df = QA)) ``` Now check also the pixel reliability information in Table 4 available via the following [link to the MODIS VI User Guide c6 version](https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf). This is important to understand how this works for the following question! @@ -198,14 +198,17 @@ plot(VI_m,1) # plot cleaned NDVI raster ``` ```{block, type="alert alert-success"} -> **Question 2**: Now what would happen if you would only work with "good" quality data? Include the new plot of the cleaned NDVI raster that only includes *good quality data" as an answer! +> **Question 2**: Now what would happen if you would only work with "good" quality data (`QA_r == 0`)? Compare the plots with both good and marginal data to the plots with only good data. ``` ```{block type="alert alert-info"} *Important*: Continue the exercise with "good and marginal data quality" as defined via the above code section in the tutorial (just rerun this section). ``` -You can (optional!) extract data from the cleaned VI raster brick via the `click` function: +You can (optional!) extract data from the cleaned VI raster brick via the `click` function. +Note: it will wait for you to click on the plot to select the pixel whose values it will print! +While it's waiting for you, no other code can be run. +If you want to cancel without clicking, press the Esc button on your keyboard. ```{r, eval=FALSE} # extract data from the cleaned raster for selected pixels click(VI_m, id=TRUE, xy=TRUE, cell=TRUE, n=1) @@ -217,7 +220,7 @@ click(VI_m, id=TRUE, xy=TRUE, cell=TRUE, n=1) ```{r, warning=FALSE, results='asis', eval=FALSE} library(leaflet) -r <- raster(VI_m,1) +r <- raster(VI_m)[[1]] # Select only the first layer (as a RasterLayer) pal <- colorNumeric(c("#ffffff", "#4dff88", "#004d1a"), values(r), na.color = "transparent") @@ -234,7 +237,7 @@ Below we extract the data from the raster as a vector and create a time series u ## the dimensions of the raster are: 33x33 px <- 78 # pixel number; adjust this number to select the center pixel -tspx <- timeser(as.vector(VI_m[px]),as.Date(names(VI_m), "X%Y.%m.%d")) # convert pixel "px" to a time series +tspx <- timeser(unlist(VI_m[px]),as.Date(names(VI_m), "X%Y.%m.%d")) # convert pixel "px" to a time series plot(tspx, main = 'NDVI') # NDVI time series cleaned using the "reliability information" ``` @@ -251,12 +254,13 @@ plot(bfm1) ``` ```{block, type="alert alert-success"} -> **Question 3**: A valid data range of NDVI in vegetation is between 0 and 1. So we should actually set NDVI values smaller than 0 to NA, as those are very likely to be invalid values. Now, do that for pixel nr. 33 and run bfastmonitor on this further cleaned NDVI time series. You can use this R code snippet for that -`tspx[tspx < 0] <- NA` and include the new bfm1 plot in your answer. Describe shortly what happens. Is this type of cleaning influencing the `bfastmonitor` analysis? Yes or No, and explain in one sentence. +> **Question 3**: A valid data range of NDVI in vegetation is between 0 and 1. So we should actually set NDVI values smaller than 0 to NA, as those are very likely to be invalid values. +Now, do that for pixel nr. 33 and run bfastmonitor on this further cleaned NDVI time series. You can use this R code snippet for that +`tspx[tspx < 0] <- NA` and check the result by plotting the time series before and after. Describe what happens. Is this type of cleaning influencing the `bfastmonitor` analysis? ``` ```{block, type="alert alert-success"} -> **Question 4**: Now check if you detect a break in the center pixel (see pixel numbering scheme). Yes or no, include the `plot(bmf1)` result for the center pixel in your answer. +> **Question 4**: Now check if you detect a break in the center pixel (see pixel numbering scheme in the website above question 1). Check the `plot(bmf1)` result for the center pixel. ``` diff --git a/index.html b/index.html index 141eb18..487d8a7 100644 --- a/index.html +++ b/index.html @@ -22,7 +22,7 @@ padding-top: 50px; padding-bottom: 50px; hyphens: auto; - word-wrap: break-word; + overflow-wrap: break-word; text-rendering: optimizeLegibility; font-kerning: normal; } @@ -31,6 +31,9 @@ font-size: 0.9em; padding: 1em; } + h1 { + font-size: 1.8em; + } } @media print { body { @@ -93,6 +96,7 @@ pre code { padding: 0; overflow: visible; + overflow-wrap: normal; } .sourceCode { background-color: transparent; @@ -134,6 +138,12 @@ #TOC li { list-style: none; } + #TOC ul { + padding-left: 1.3em; + } + #TOC > ul { + padding-left: 0; + } #TOC a:not(:hover) { text-decoration: none; } @@ -1003,7 +1013,7 @@
@@ -1023,47 +1033,106 @@

14 June, 2021

} -

MODIS based time series analysis using BFAST Monitor and BFAST Lite

-

WUR logo

-

Jan Verbesselt, Dainius Masiliūnas

+

MODIS +based time series analysis using BFAST Monitor and BFAST Lite

+

WUR logo

+

Jan Verbesselt, Dainius +Masiliūnas

-

14-06-2021

-
-

This document explains how to use the R scripting language for downloading MODIS data and analyzing its time series within R. By the end of the tutorial, you will be able to download and preprocess MODIS data, and apply a time series break detection algorithm to automatically determine changes in time series. For this time series analysis demonstration it is not required to know R details, we only use R for some practical demonstration of its great potential. In this exercise we will automatically download MODIS data for specific locations around the world. You can then apply, at your own choice, one of the two algorithms for break detection: BFAST Monitor or BFAST Lite.

-

Introduction to MODIS data and online analysis

-

The MODIS satellite data that we will download for this time series analysis exercise is available from the MODIS Subsetting Website. MODIS data is made available for specific locations.

-

More specifically, we will look at the MODIS product called MOD13Q1 which are global 16-day images at a spatial resolution of 250 m. Each image contains several bands; i.e. blue, red, and near-infrared reflectances, centered at 469 nanometers, 645 nanometers, and 858 nanometers, respectively. These bands are used to determine the MODIS vegetation indices.

-

The MODIS Normalized Difference Vegetation Index (NDVI) complements NOAA’s Advanced Very High Resolution Radiometer (AVHRR) NDVI products and provides continuity for applications requiring long time series. MODIS also includes a new Enhanced Vegetation Index (EVI) that minimizes canopy background variations and maintains sensitivity over dense vegetation conditions. The EVI also uses the blue band to remove residual atmosphere contamination caused by smoke, sub-pixel and thin clouds. The MODIS NDVI and EVI products are computed from atmospherically corrected surface reflectances that have been masked for water, clouds, heavy aerosols, and cloud shadows.

-

Vegetation indices are used for global monitoring of vegetation conditions and are used in products displaying land cover and its changes. We will work with the MODIS NDVI band of the MOD13Q1 product. More information about this MODIS data set can be found via the MODIS product table.

-

Now go to the MODIS Land Products:

+

18-03-2022

+
+

This document explains how to use the R scripting language for +downloading MODIS data and analyzing its time series within R. By the +end of the tutorial, you will be able to download and preprocess MODIS +data, and apply a time series break detection algorithm to automatically +determine changes in time series. For this time series analysis +demonstration it is not required to know R details, we only use R for +some practical demonstration of its great potential. In this exercise we +will automatically download MODIS data for specific locations around the +world. You can then apply, at your own choice, one of the two algorithms +for break detection: BFAST Monitor or BFAST Lite.

+

Introduction to +MODIS data and online analysis

+

The MODIS satellite data that we will download for this time series +analysis exercise is available from the MODIS Subsetting Website. MODIS +data is made available for specific locations.

+

More specifically, we will look at the MODIS product called MOD13Q1 +which are global 16-day images at a spatial resolution of 250 +m. Each image contains several bands; i.e. blue, red, and +near-infrared reflectances, centered at 469 nanometers, 645 nanometers, +and 858 nanometers, respectively. These bands are used to determine the +MODIS vegetation indices.

+

The MODIS Normalized Difference Vegetation Index (NDVI) complements +NOAA’s Advanced Very High Resolution Radiometer (AVHRR) NDVI products +and provides continuity for applications requiring long time series. +MODIS also includes a new Enhanced Vegetation Index (EVI) that minimizes +canopy background variations and maintains sensitivity over dense +vegetation conditions. The EVI also uses the blue band to remove +residual atmosphere contamination caused by smoke, sub-pixel and thin +clouds. The MODIS NDVI and EVI products are computed from +atmospherically corrected surface reflectances that have been masked for +water, clouds, heavy aerosols, and cloud shadows.

+

Vegetation indices are used for global monitoring of vegetation +conditions and are used in products displaying land cover and its +changes. We will work with the MODIS NDVI band of the MOD13Q1 product. +More information about this MODIS data set can be found via the MODIS product +table.

+

Now go to the MODIS Land +Products:

-Question 1: What are the three main land cover types close to the center of the site? (See the integrated legend in the maps below.) Which land cover types are correct and which are not? +Question 1: What are the three main land cover types +close to the center of the site? (See the integrated legend in the maps +below.) Which land cover types are correct and which are not?

-

Preprocessing: automated MODIS NDVI download and analysis via R

-

We will download the MODIS data for the Loobos Site via R and process the data for one location to detect changes within the time series.

+

Preprocessing: +automated MODIS NDVI download and analysis via R

+

We will download the MODIS data for the Loobos Site via R and process +the data for one location to detect changes within the time series.

-

Getting started: install packages and load the necessary functions for MODIS data analysis

-

Now we are ready to get started with the MODIS time series analysis exercise in R!

-First, choose your working directory (i.e. a folder on your hard drive) to which the MODIS data will be downloaded and where you will save your R script. +

Getting +started: install packages and load the necessary functions for MODIS +data analysis

+

Now we are ready to get started with the MODIS time series analysis +exercise in R!

+First, choose your working directory (i.e. a folder on your hard drive) +to which the MODIS data will be downloaded and where you will save your +R script.

-Protip: Do not hard code setwd() in your R script as the path might be different on another computer, and therefore your script will not be fully reproducible by others. +Protip: Do not hard code setwd() in your R +script as the path might be different on another computer, and therefore +your script will not be fully reproducible by others.

-

In RStudio you can set your working directory this way, if you have saved your R script to your working directory:

+

In RStudio you can set your working directory this way, if you have +saved your R script to your working directory:

Check your working directory by

-

The necessary add-on packages need to be installed within R before loading the package using the library() function. Below we define a helper function that installs the R package if it is not installed yet, and then also loads it using the library function:

+

The necessary add-on packages need to be installed within R before +loading the package using the library() function. Below we +define a helper function that installs the R package if it is not +installed yet, and then also loads it using the library +function:

-

Loading extra function timeser() to create a time series object in R:

+

Loading extra function timeser() to create a time series +object in R:

@@ -1197,7 +1288,9 @@

-Question 2: Now what would happen if you would only work with “good” quality data? Include the new plot of the cleaned NDVI raster that only includes *good quality data" as an answer! +Question 2: Now what would happen if you would only +work with “good” quality data? Include the new plot of the cleaned NDVI +raster that only includes *good quality data” as an answer!

@@ -1205,11 +1298,14 @@

-Important: Continue the exercise with “good and marginal data quality” as defined via the above code section in the tutorial (just rerun this section). +Important: Continue the exercise with “good and marginal data +quality” as defined via the above code section in the tutorial (just +rerun this section).

-You can (optional!) extract data from the cleaned VI raster brick via the click function: +You can (optional!) extract data from the cleaned VI raster brick via +the click function:
@@ -1239,7 +1336,8 @@

-Below we extract the data from the raster as a vector and create a time series using the timeser function: +Below we extract the data from the raster as a vector and create a time +series using the timeser function:
-

Now we are ready to detect breaks in the time series! You can now choose: either use BFAST Monitor (“Option 1”, the following section) to detect a single break at the end of the time series, or use BFAST Lite (“Option 2”, the section after that) to detect all breaks in the middle of the time series. If you are interested, you can do both, but it’s not necessary to answer both sets of questions.

-

Option 1: detect break at the end of the time series with BFAST Monitor

-Now we apply the bfastmonitor function using a trend + harmon model with order 3 for the harmonics (i.e. seasonality modelling): + + + + +

Now we are ready to detect breaks in the time series! You can now +choose: either use BFAST Monitor (“Option 1”, the following section) to +detect a single break at the end of the time series, or use BFAST Lite +(“Option 2”, the section after that) to detect all breaks in the middle +of the time series. If you are interested, you can do both, but it’s not +necessary to answer both sets of questions.

+

Option +1: detect break at the end of the time series with BFAST Monitor

+Now we apply the bfastmonitor function using a +trend + harmon model with order 3 for the +harmonics (i.e. seasonality modelling):
@@ -1275,7 +1384,15 @@

-Question 3: A valid data range of NDVI in vegetation is between 0 and 1. So we should actually set NDVI values smaller than 0 to NA, as those are very likely to be invalid values. Now, do that for pixel nr. 33 and run bfastmonitor on this further cleaned NDVI time series. You can use this R code snippet for that tspx[tspx < 0] <- NA and include the new bfm1 plot in your answer. Describe shortly what happens. Is this type of cleaning influencing the bfastmonitor analysis? Yes or No, and explain in one sentence. +Question 3: A valid data range of NDVI in vegetation is +between 0 and 1. So we should actually set NDVI values smaller than 0 to +NA, as those are very likely to be invalid values. Now, do that for +pixel nr. 33 and run bfastmonitor on this further cleaned NDVI time +series. You can use this R code snippet for that tspx[tspx < 0] +<- NA and include the new bfm1 plot in your answer. Describe +shortly what happens. Is this type of cleaning influencing the +bfastmonitor analysis? Yes or No, and explain in one +sentence.

@@ -1284,7 +1401,9 @@

-Question 4: Now check if you detect a break in the center pixel (see pixel numbering scheme). Yes or no, include the plot(bmf1) result for the center pixel in your answer. +Question 4: Now check if you detect a break in the +center pixel (see pixel numbering scheme). Yes or no, include the +plot(bmf1) result for the center pixel in your answer.

@@ -1293,12 +1412,17 @@

-Question 5: Now check pixel 78. What happens if you use a different formula in bfastmonitor for the pixel? For example, response ∼ trend. Explain what happens and also in which sense is 2019 an abnormal year? Include the bfastmonitor plots in your answer. +Question 5: Now check pixel 78. What happens if you use +a different formula in bfastmonitor for the pixel? +For example, response ∼ trend. Explain what happens and +also in which sense is 2019 an abnormal year? Include the bfastmonitor +plots in your answer.

-

Let’s run the bfastmonitor code on the full raster brick spatially using the calc function:

+

Let’s run the bfastmonitor code on the full raster +brick spatially using the calc function:

@@ -1327,7 +1452,9 @@

-Question 6: Explain what you see in these two plots (created with the above code section). Do you think these plots imply major changes for the loobos site? Justify why. +Question 6: Explain what you see in these two plots +(created with the above code section). Do you think these plots imply +major changes for the loobos site? Justify why.

@@ -1336,12 +1463,21 @@

-Question 7: Now try to detect a change in 2019 for another site in the world available via the https://modis.ornl.gov/sites/ website. In your answer, mention (a) site, (b) bfastmonitor settings, (c) landcover type of that pixel and the pixel number, (d) include the the time.of.break and magnitude.of.change plots in your answer, (e) briefly discuss the main differences and/or similarities (in magnitude) of your site with the loobos site. +Question 7: Now try to detect a change in 2019 for +another site in the world available via the https://modis.ornl.gov/sites/ +website. In your answer, mention (a) site, (b) bfastmonitor settings, +(c) landcover type of that pixel and the pixel number, (d) include the +the time.of.break and magnitude.of.change plots in your answer, (e) +briefly discuss the main differences and/or similarities (in magnitude) +of your site with the loobos site.

-Here is example R code to get the pixel number. When you run click(), click in the plot on a pixel with a break (i.e. where an estimated time of break is available). +Here is example R code to get the pixel number. When you run +click(), click in the plot on a pixel with a break +(i.e. where an estimated time of break is available).
-Here we selected one pixel, and do the bfastmonitor analysis for that pixel: +Here we selected one pixel, and do the bfastmonitor +analysis for that pixel:
-

Option 2: detecting all breaks in the middle of a time series with BFAST Lite

-

BFAST Monitor is made to detect the first break at the end of a time series. If you need to detect more than one break, then you need to use a different algorithm: either BFAST or BFAST Lite. In this example, we will run BFAST Lite (a lightweight version of BFAST that can handle NA values) on the same data as we did above.

-

Let’s run the function bfastlite() on our data from the previous steps. Setting the parameter breaks to BIC results in more liberal detection of breaks compared to the default, LWZ.

+ + + + +

Option +2: detecting all breaks in the middle of a time series with BFAST +Lite

+

BFAST Monitor is made to detect the first break at the end of a +time series. If you need to detect more than one break, then you +need to use a different algorithm: either BFAST or BFAST Lite. In this +example, we will run BFAST Lite (a lightweight version of BFAST that can +handle NA values) on the same data as we did above.

+

Let’s run the function bfastlite() on our data from the +previous steps. Setting the parameter breaks to +BIC results in more liberal detection of breaks compared to +the default, LWZ.

-

Why did the model decide to place one break? Let’s take a look at the statistics:

+

Why did the model decide to place one break? Let’s take a look at the +statistics:

-

The residual sum of squares keeps decreasing with more breakpoints, as that makes the model more flexible and thus fit the data better, but BIC and other information criteria apply a penalty for adding more breaks than really necessary, thus limiting the number of false detections. In this case, both BIC and LWZ agree that 1 is the optimal number of breaks.

-

We can visually see that the identified break is indeed quite significant. We can also look at it from a statistical point of view by using the function magnitude():

+

The residual sum of squares keeps decreasing with more breakpoints, +as that makes the model more flexible and thus fit the data better, but +BIC and other information criteria apply a penalty for adding more +breaks than really necessary, thus limiting the number of false +detections. In this case, both BIC and LWZ agree that 1 is the optimal +number of breaks.

+

We can visually see that the identified break is indeed quite +significant. We can also look at it from a statistical point of view by +using the function magnitude():

-

This shows that the difference between the actual observations and the predictions of the models on both sides of the break, when extrapolated to the other side, is fairly high (0.22 RMSD and MAD, thus at the magnitude of 0.22 NDVI units).

+

This shows that the difference between the actual observations and +the predictions of the models on both sides of the break, when +extrapolated to the other side, is fairly high (0.22 RMSD and MAD, thus +at the magnitude of 0.22 NDVI units).

-Question 3: Try changing the parameter breaks in the bfastlite() function to LWZ and to an integer number, and rerun the code. What is the result? How is it different from what we had above? Explain in one sentence. +Question 3: Try changing the parameter +breaks in the bfastlite() function to +LWZ and to an integer number, and rerun the code. What is +the result? How is it different from what we had above? Explain in one +sentence.

@@ -1556,7 +1725,9 @@

-Question 4: Now check if you detect a break in the center pixel (see pixel numbering scheme). Yes or no, include the plot(breaks) result for the center pixel in your answer. +Question 4: Now check if you detect a break in the +center pixel (see pixel numbering scheme). Yes or no, include the +plot(breaks) result for the center pixel in your answer.

@@ -1565,12 +1736,19 @@

-Question 5: Now check pixel 82 again. What happens if you use a different formula in bfastlite for the pixel? For example, response ∼ trend. Explain what happens and also in which sense is 2016 an abnormal year? Include the bfastlite plots in your answer. +Question 5: Now check pixel 82 again. What happens if +you use a different formula in bfastlite for the +pixel? For example, response ∼ trend. Explain what happens +and also in which sense is 2016 an abnormal year? Include the bfastlite +plots in your answer.

-

Let’s try to run BFAST Lite on the whole raster now. There is one problem with that: since the number of breaks is variable, we don’t have a variable number of layers. Thus, let’s plot only the break that is the most significant, so that we get two layers as output.

+

Let’s try to run BFAST Lite on the whole raster now. There is one +problem with that: since the number of breaks is variable, we don’t have +a variable number of layers. Thus, let’s plot only the break that is the +most significant, so that we get two layers as output.

-

Plot the results (optionally compare with BFAST Monitor results above):

+

Plot the results (optionally compare with BFAST Monitor results +above):

@@ -1623,7 +1803,9 @@

-Question 6: Explain what you see in these two plots (created with the above code section). Do you think these plots imply major changes for the loobos site? Justify why. +Question 6: Explain what you see in these two plots +(created with the above code section). Do you think these plots imply +major changes for the loobos site? Justify why.

@@ -1632,15 +1814,27 @@

-Question 7: Now try to detect a change for another site in the world available via the https://modis.ornl.gov/sites/ website. In your answer, mention (a) site, (b) bfastlite settings, (c) landcover type of that pixel and the pixel number, (d) include the the time.of.break and magnitude.of.change plots in your answer, (e) briefly discuss the main differences and/or similarities (in magnitude) of your site with the loobos site. +Question 7: Now try to detect a change for another site +in the world available via the https://modis.ornl.gov/sites/ +website. In your answer, mention (a) site, (b) bfastlite settings, (c) +landcover type of that pixel and the pixel number, (d) include the the +time.of.break and magnitude.of.change plots in your answer, (e) briefly +discuss the main differences and/or similarities (in magnitude) of your +site with the loobos site.

More information

-

More information can be found on the BFAST website and in the BFAST papers mentioned on the website.

-

The following section gives extra information about the concept of seasonality monitoring using harmonic analysis. There are no questions from this section, it’s for your interest only,

-

Seasonality monitoring using harmonics

+

More information can be found on the BFAST website and in the BFAST +papers mentioned on the website.

+

The following section gives extra information about the concept of +seasonality monitoring using harmonic analysis. There are no questions +from this section, it’s for your interest only,

+

Seasonality monitoring +using harmonics

-





Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

+





+Creative Commons License
This +work is licensed under a +Creative +Commons Attribution-ShareAlike 4.0 International License.

@@ -1693,7 +1893,8 @@

Seasonality monitoring using har