MODISTools2/inst/doc/UsingMODISTools.Rnw

\documentclass[11pt]{article}
\usepackage[top=2cm, bottom=3cm, left=2cm, right=2cm]{geometry}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}            % /eqref
\usepackage{url}
\usepackage{hyperref}
\usepackage[none]{hyphenat}     % No hyphens
\usepackage[noae]{Sweave}
%\VignetteIndexEntry{Using MODISTools}
\newcommand{\code}[1]{\texttt{#1}}

\begin{document}

\sloppy    % Prevent hyphenated words running into margins
\SweaveOpts{concordance=TRUE, width=6, height=6}
\setkeys{Gin}{width=0.5\textwidth}
\title{Using MODISTools
       (\Sexpr{packageDescription('MODISTools', fields='Version')})}
\author{Sean Tuck}
\date{\Sexpr{packageDescription('MODISTools', fields='Date')}}

\maketitle
\tableofcontents

<<echo=FALSE>>=
library(MODISTools)

# Makes copy-paste much less painful
options(continue = ' ')
options(width = 90)
options(prompt = '> ')

options(SweaveHooks = list(fig=function() par(mgp=c(2.5,1,0),
                                              mar=c(4,4,2,1),
                                              oma=c(0,0,1,0),
                                              cex.main=0.8)))
@

\section{Introduction}
The MODISTools R package is a set of tools for downloading and working with NASA's MODIS remotely-sensed data. The package retrieves data from the LP DAAC data archive, via their SOAP web service. Functions download data as a batch process, and save subsets in text files that can be returned to at a later date. Additional functions can provide summaries of this data and prepare the data to a format ready for application in R; if you have other data that you wish to relate MODIS data to, downloaded data can be appended to your original dataset. Other ancillary functions can help to get input arguments into the correct format.

This vignette provides a worked example for using MODISTools. A dataset of time-series -- lat-long coordinates with start and end dates -- to collect MODIS data for, will be used to show a complete workflow for how someone might use MODISTools. We will prepare input information for a subset request, download subsets of Enhanced Vegetation Index (EVI) and land cover data for the specified locations, and process these data to analyse land processes at these locations. Note that you will need an internet connection to run this worked example yourself, and that it will download files to your computer.

\section{Format the data}
We have some coordinates that we would like to extract MODIS data for. But the coordinates are not in the correct format. We need to make sure the coordinates we input for our subset request are in the WGS-1984 coordinate system, and are in decimal degrees format.
<<>>=
data(ConvertExample)
ConvertExample
@
These coordinates are WGS-1984 coordinates, but they are not in decimal degrees. We can use \code{ConvertToDD} to fix this.
<<>>=
modis.subset <-
  ConvertToDD(XY = ConvertExample, LatColName = "lat", LongColName = "long")
modis.subset <- data.frame(lat = modis.subset[ ,1], long = modis.subset[ ,2])
modis.subset
@
What we also need to retrieve a time-series of MODIS data for these locations are dates. End dates for the time-series, and preferably start dates too. If we don't have start dates we can ask for a set number of years for each location instead. Let's retrieve data between 2003 and 2006. The dates can be specified as years or in POSIXlt date-time class (see \code{?POSIXlt}). In this case we can just use years.
<<>>=
modis.subset$start.date <- rep(2003, nrow(modis.subset))
modis.subset$end.date <- rep(2006, nrow(modis.subset))
@
That's all we need! Let's download our EVI data first.

\section{Download the data}
\subsection{Specifying a subset request}
The shortname code for the EVI product is ``MOD13Q1". We can check the codes for all the products available using \code{GetProducts}, and we can find the shortname codes for all data bands within each product using \code{GetBands}.
<<eval=FALSE>>=
GetProducts()
@
<<echo=FALSE>>=
c("MCD12Q1", "MCD12Q2", "MCD43A1", "MCD43A2", "MCD43A4", "MOD09A1",
"MOD11A2", "MOD13Q1", "MOD15A2", "MOD15A2GFS", "MOD16A2", "MOD17A2_51",
"MOD17A3", "MYD09A1", "MYD11A2", "MYD13Q1", "MYD15A2")
@
<<eval=FALSE>>=
GetBands(Product = "MOD13Q1")
@
<<echo=FALSE>>=
c("250m_16_days_blue_reflectance", "250m_16_days_MIR_reflectance",
"250m_16_days_NIR_reflectance", "250m_16_days_pixel_reliability",
"250m_16_days_red_reflectance", "250m_16_days_relative_azimuth_angle",
"250m_16_days_sun_zenith_angle", "250m_16_days_view_zenith_angle",
"250m_16_days_VI_Quality", "250m_16_days_NDVI",
"250m_16_days_EVI", "250m_16_days_composite_day_of_the_year")
@
We will download EVI data at 250m pixel resolution, which is available at 16-day intervals. The shortname code for this data band is 250m\_16\_days\_EVI. We will collect quality control data for these pixels too, which is available from the 250m\_16\_days\_pixel\_reliability band (and 250m\_16\_days\_VI\_Quality too).

We can check that the time-series of MODIS data we want is available for this data product by retrieving the dates for all available time-steps.
<<eval=FALSE>>=
GetDates(Product = "MOD13Q1", Lat = modis.subset$lat[1], Long = modis.subset$long[1])
@
The time-period available for the Vegetation Indices product covers 2003-2006 (the maximum shown is at the time this vignette was built), so we can proceed. When we download we also need to decide how large we want the tiles of data for each location to be. We specify this by entering the distance (km) above and below in each direction away from the central pixel, where the input coordinate is located, and then doing the same for left and right. The input must be whole km (integers) for each direction. As an example, if we specify \code{Size=c(1,1)} for this EVI data at 250m pixel resolution, it will retrieve a 9x9 pixel tile for each location, centred on the input coordinate. The tiles this size will be downloaded at the locations for each time-step that falls between the start and end dates. \code{Size=c(0,0)} would specify only the central pixel. The maximum size tile surrounding a location is \code{Size=c(100,100)}.

\subsection{MODISSubsets}
The download will write the MODIS data to ASCII files for each location subset specified. We can specify the directory that we would like to save downloaded files in, using the \code{SaveDir} argument below. In the code below, downloaded files will be written to your working directory; if you would prefer the files to be written elsewhere change \code{SaveDir}. But we will access these files later, so remember to request the files from the same directory.
<<eval=FALSE>>=
MODISSubsets(LoadDat = modis.subset, Products = "MOD13Q1",
             Bands = c("250m_16_days_EVI", "250m_16_days_pixel_reliability"),
             Size = c(1,1))
@
Each ASCII file is a different subset location. In each ASCII file, each row is a different time-step in the time-series. If multiple data bands have been downloaded for this subset, they will all be contained in the same ASCII file for that subset.

Here is an example of the strings of data that are downloaded for pixels at each time-step and data band:
<<eval=FALSE>>=
subset.string <- read.csv(list.files(pattern = ".asc")[1],
                                header = FALSE, as.is = TRUE)
subset.string[1, ]
@
<<echo=FALSE>>=
subset.string <- read.csv(paste("./MODISSubsetsMOD13Q1/",
                                list.files(path = "./MODISSubsetsMOD13Q1", pattern = ".asc")[1]
                                , sep = ""), header = FALSE, as.is = TRUE)
subset.string[1, ]
@
A download log file will also be written, displaying all the unique subsets found in the dataset, and confirmation of download success for each.

\subsection{MODISTransects}
Alternatively, we may want transects of MODIS data. This is easily done by specifying start and end points for transects, with unique IDs for each transect and then calling \code{MODISTransects} instead of using \code{MODISSubsets}. To try this out, use the example given in the \code{MODISTransects} help documentation.

\section{Process the data}
\subsection{MODISSummaries}
Now we have downloaded the EVI data, we can find average each pixel over time, to produce one tile of mean EVI pixels at each subset location. We can use \code{MODISSummaries} for this. The function will also take this processed data and append it to your original files containing all the subset information (\code{modis.subset}). This will write two files to the specified directory. We downloaded quality control data for each pixel alongside our EVI data, so \code{MODISSummaries} can also check for poor quality and missing data. These data will be removed and replaced with \code{NA}s. The threshold for defining what is good and poor quality is set by the user: the scores for highest quality is 0, and the score for lowest quality is 3 or 5, depending on the data band. To see how quality control information is defined for each data type, go to the \href{https://daacweb-dev.ornl.gov/MODIS/MODIS-menu/products.html}{\bf{MODIS Products Table}}. We need to specify the range of valid data for EVI, the value that denotes missing data, and the scale factor that is applied to the data, which are all available from the same web page.
<<eval=FALSE>>=
MODISSummaries(LoadDat = modis.subset, Product = "MOD13Q1", Bands = "250m_16_days_EVI",
               ValidRange = c(-2000,10000), NoDataFill = -3000, ScaleFactor = 0.0001,
               QualityScreen = TRUE, QualityBand = "250m_16_days_pixel_reliability",
               QualityThreshold = 0)
@
If you want to screen data for quality without all the other things that \code{MODISSummaries} does, you can call the more general \code{QualityCheck}, which is an internal function for \code{MODISSummaries}.

\subsection{ExtractTile}
Also, if large subset tiles are downloaded for each location, there may be times when we want to extract a smaller tile from within this subset, rather than downloading again to retrieve the nested data we want. This can be done using \code{ExtractTile}. We will use the file just written from our call to \code{MODISSummaries}, retrieve the smaller subset we want, and arrange them into tiles to compare before and after.
<<eval=FALSE>>=
TileExample <- read.csv(list.files(pattern = "MODIS_Data"))
TileExample <- TileExample[ ,which(grepl("250m_16_days_EVI", names(TileExample)))]
@
<<echo=FALSE>>=
TileExample <- read.csv(paste("./MODISSummaries/",
                              list.files(path = "./MODISSummaries/",
                                         pattern = "Data"), sep = ""))
TileExample <- TileExample[ ,which(grepl("250m_16_days_EVI", names(TileExample)))]
@
Pixels in a tile are on the same row. See that using \code{ExtractTile} takes away some of the columns.
<<>>=
dim(TileExample)
dim(ExtractTile(Data = TileExample, Rows = c(9,2), Cols = c(9,2), Grid = FALSE))
head(ExtractTile(Data = TileExample, Rows = c(9,2), Cols = c(9,2), Grid = FALSE),
     n = 2)
@
We can look at the first subset and arrange the pixels into a tile to visually show what \code{ExtractTile} has done.
<<>>=
matrix(TileExample[1, ], nrow = 9, ncol = 9, byrow = TRUE)
ExtractTile(Data = TileExample, Rows = c(9,2), Cols = c(9,2), Grid = TRUE)[ , ,1]
@
Arrangement of the pixels into tiles this way can be optionally set with a call to \code{ExtractTile}. The order for the strings of pixel data in the downloaded ASCII files is by row, so \code{matrix(..., byrow=TRUE)} can arrange the pixels correctly (see above).

\subsection{LandCover}
Let's do the same as above but download data on land cover classes for the same subsets.
<<eval=FALSE>>=
dir.create('./LandCover')
setwd('./LandCover')
MODISSubsets(LoadDat = modis.subset, Product = "MCD12Q1", Bands = "Land_Cover_Type_1",
             Size = c(1,1))
@
We can use \code{LandCover} to retrieve some summaries of land cover in each tile. This will tell us the most common land cover type, the total number of distinct land cover types, and Simpson's D and evenness measures to express landscape diversity and heterogeneity in these tiles. Let's retrieve these summaries from the land cover subset files we just downloaded.
<<eval=FALSE>>=
LandCover(Band = "Land_Cover_Type_1")

land.summary <- read.csv(list.files(pattern = "MODIS_Land_Cover_Summary"))
head(land.summary)
@
<<echo=FALSE>>=
land.summary <- read.csv(paste("./LandCover/",
                               list.files(path = "./LandCover/",
                                          pattern = "LandCoverSummary"),
                               sep = ""))
head(land.summary)
@

\end{document}