land_cover.Rmd

---
title: "Land Cover"
output:
  html_document:
    toc: true
    toc_float: true
    number_sections: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```

#### Data set details

|||
| ----------- | ----------- |
| **Data set description:** | Land cover classification in Zimbabwe |
| **Source:**   | [Moderate Resolution Imaging Spectroradiometer (MODIS)](https://modis.gsfc.nasa.gov/data/dataprod/) |
| **Details on the retrieved data:**   | Land cover classification ([Land Cover Products](https://modis.gsfc.nasa.gov/data/dataprod/mod12.php)) in Zimbabwe in 2019. |
| **Spatial and temporal resolution:**  | Annually averaged data at a 500 meter spatial resolution. |


## MODIS (Moderate Resolution Imaging Spectroradiometer)

MODIS is an instrument aboard the Terra and Aqua satellites, which orbits the entire Earth every 1-2 days, acquiring data at different spatial resolutions. The data acquired by MODIS describes features of the land, oceans and the atmosphere. A complete list of MODIS data products can be found on the [MODIS website](https://modis.gsfc.nasa.gov/data/dataprod/).

## Downloading Land Cover data using `MODIStsp`

`MODIStsp` is an R package for downloading and preprocessing time series of raster data from MODIS data products.
The package's name is an acronym for 'MODIS Time Series Processing'.
This tutorial will focus on downloading and visualising land cover data, but the same process can be followed with other MODIS data products as well.

### Installing `MODIStsp`

The `MODIStsp` package can be downloaded from CRAN as follows.

```{r install-MODIStsp, eval=FALSE}
install.packages("MODIStsp")
```

### Identifying the MODIS data product

The first step of downloading data is to identify which [MODIS data product](https://modis.gsfc.nasa.gov/data/dataprod/) to use.

This tutorial will use the [Land Cover Products](https://modis.gsfc.nasa.gov/data/dataprod/mod12.php), which uses different classification schemes to identify types of land cover. The primary land cover scheme identifies 17 classes defined by the **International Geosphere-Biosphere Programme (IGBP)** as detailed in their [user guide](https://lpdaac.usgs.gov/documents/101/MCD12_User_Guide_V6.pdf).This product contains annually averaged data, collected bi-annually allowing for hemispheric differences in the growing seasons, and enabling the product to capture two growth cycles if necessary.

The product IDs for each of these products can also be found on the [data product page](https://modis.gsfc.nasa.gov/data/dataprod/mod13.php). 

This tutorial will use the 'Land Cover Type Yearly L3 Global 500m' product with the product IDs MCD12Q1 (Combined Aqua and Terra Product ID)

### Retrieving MODIS layers for a product

The product layers (original MODIS layers, quality layers and spectral indexes) available for a given product can be retrieved using the following function.

```{r modis-get-productlayers, attr.output='style="max-height: 200px;"'}
library(MODIStsp)
MODIStsp_get_prodlayers("MCD12Q1")
```

Note how the `\$bandfullnames` define each of the `\$bandnames`.

### The `MODIStsp()` function

`MODIStsp()` is the main function of the `MODIStsp` package, and allows us to download MODIS data products. While this is a very comprehensive function and we only a very few of its arguments in this tutorial, the entire list of arguments can be found in the [MODIStsp documentation](https://docs.ropensci.org/MODIStsp/reference/MODIStsp.html).

The `MODIStsp()` function provides two ways of downloading data, namely, through a GUI (interactive) or through an R script (non-interactive). This tutorial will focus on the non-interactive execution.

### Downloading IGBP land cover classification data

To download the IGBP land cover classification data in Zimbabwe, first we download the boundary of Zimbabwe with the `geoboundaries()` function from the `rgeoboundaries` package and save it on our computer.

```{r downloading-shapefile, eval=FALSE}
# remotes::install_github("wmgeolab/rgeoboundaries")
# install.packages("sf")
library(rgeoboundaries)
library(sf)
# Downloading the country boundary of Zimbabwe
map_boundary <- geoboundaries("Zimbabwe")
# Defining filepath to save downloaded spatial file
spatial_filepath <- "LandCoverData/zimbabwe.shp"
# Saving downloaded spatial file on to our computer
st_write(map_boundary, paste0(spatial_filepath))
```

Then we use the `MODIStsp()` function to download the IGBP land cover classification data.

To download data in Zimbabwe, we use the boundary of Zimbabwe we downloaded. So in the `MODIStsp()` function we set the `spatmeth` argument as "file" and set the `spafile` argument as the path of the map we saved.

Since IGBP land cover classification data is provided annually and dated on the 1st of January each year, the `start-date` and `end_date` arguments should contain the 1st of January of each year for which we want the data to be downloaded.

The `start-date` and `end_date` arguments define the period for which we want the data to be downloaded.

Note that this tutorial uses a test username and password. The `user` and `password` arguments should be the username and password corresponding to your [earthdata](https://urs.earthdata.nasa.gov/) credentials.

```{r downloading-data, eval=FALSE}
library(MODIStsp)
MODIStsp(gui             = FALSE,
         out_folder      = "LandCoverData",
         out_folder_mod  = "LandCoverData",
         selprod         = "LandCover_Type_Yearly_500m (MCD12Q1)",
         bandsel         = "LC1", 
         user            = "mstp_test" ,
         password        = "MSTP_test_01",
         start_date      = "2019.01.01", 
         end_date        = "2019.12.31", 
         verbose         = FALSE,
         spatmeth        = "file",
         spafile         = spatial_filepath,
         out_format      = "GTiff")
```

## Understanding the downloaded data

The downloaded files are saved in subfolders within the defined output folder.

A separate subfolder is created for each processed original MODIS layer, Quality Indicator or Spectral Index with an image for each processed date. The images will be placed in the following folder structure and named using the following naming convention.

`<defined_out_folder>/<shape_file_name>/<product_name>/<layer_name>/<prodcode>_<layername>_<YYYY>_<day_of_year>.<extension>`

## Visualising land cover in Zimbabwe

The following example uses the `geom_raster()` function from the `ggplot2` package to visualise the downloaded IGBP land cover classification in Zimbabwe

```{r visualise-data}
# install.packages(c("sf", "raster", "here", "ggplot2", "viridis"))
# remotes::install_github("wmgeolab/rgeoboundaries")
library(rgeoboundaries)
library(sf)
library(raster)
library(here)
library(ggplot2)
library(viridis)
library(dplyr)

# Downloading the boundary of Zimbabwe
map_boundary <- geoboundaries("Zimbabwe")

# Reading in the downloaded landcover raster data
IGBP_raster <- raster(here::here("LandCoverData/zimbabwe/LandCover_Type_Yearly_500m_v6/LC1/MCD12Q1_LC1_2019_001.tif"))

# Transforming data
IGBP_raster <- projectRaster(IGBP_raster, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Cropping data
IGBP_raster <- raster::mask(IGBP_raster, as_Spatial(map_boundary))

# Converting the raster object into a dataframe and converting the IGBP classification into a factor
IGBP_df <- as.data.frame(IGBP_raster, xy = TRUE, na.rm = TRUE) %>%
  mutate(MCD12Q1_LC1_2019_001 = as.factor(round(MCD12Q1_LC1_2019_001)))
rownames(IGBP_df) <- c()
# Renaming IGBP classification levels
levels(IGBP_df$MCD12Q1_LC1_2019_001) <- c( "Evergreen needleleaf forests",
                                           "Evergreen broadleaf forests",
                                           "Deciduous needleleaf forests",
                                           "Deciduous broadleaf forests",
                                           "Mixed forests",
                                           "Closed shrublands",
                                           "Open shrublands",
                                           "Woody savannas",
                                           "Savannas",
                                           "Grasslands",
                                           "Permanent wetlands",
                                           "Croplands",
                                           "Urban and built-up lands",
                                           "Cropland/natural vegetation mosaics",
                                           "Snow and ice",
                                           "Barren",
                                           "Water bodies")
# Visualising using ggplot2
ggplot() + 
  geom_raster(data = IGBP_df,
              aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
  geom_sf(data = map_boundary, inherit.aes = FALSE, fill = NA) +
  scale_fill_viridis(name = "Land Cover Type", discrete = TRUE) +
  labs(title = "Land Cover classification in Zimbabwe",
       subtitle = "01-01-2019 - 31-12-2019",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal()
```


## References

- `MODIS` website: https://modis.gsfc.nasa.gov/

- `MODIStsp` vignette: https://cran.r-project.org/web/packages/MODIStsp/vignettes/MODIStsp.html

- `MODIStsp land cover` data product page: https://modis.gsfc.nasa.gov/data/dataprod/mod12.php


---

<p style ="color:grey; font-size: 13px;">  
Last updated: `r Sys.Date()`  
Source code: https://github.com/rspatialdata/rspatialdata.github.io/blob/main/land_cover.Rmd 
<details>
<summary><span style ="color:grey; font-size: 13px;">Tutorial was complied using: (click to expand)</span></summary>
```{r echo=FALSE}
sessionInfo()
```
</details>
</p>