---
title: "Satellite Remote Sensing"
---
```{r, echo=FALSE, message=FALSE, results='hide', purl=FALSE}
## This chunk automatically generates a text .R version of this script when running within knitr. You do not need to run this...
input = knitr::current_input() # filename of input document
output = paste(tools::file_path_sans_ext(input), 'R', sep = '.')
knitr::purl(input,output,documentation=2,quiet=T)
source("knitr_header.R")
knitr::opts_chunk$set(eval=T)
```
[ The R Script associated with this page is available here](`r output`). Download this file and open it (or copy-paste into a new script) with RStudio so you can follow along.
### Libraries
```{r,results='hide',message=FALSE}
library(raster)
library(rasterVis)
library(rgdal)
library(sp)
library(ggplot2)
library(ggmap)
library(dplyr)
library(reshape2)
library(knitr)
library(tidyr)
# New Packages
library(MODISTools)
library(gdalUtils)
library(rts)
```
## Specify directory to store data (absolute or relative to current working directory).
```{r}
download.file("http://adamwilson.us/RDataScience/09_data.zip",
destfile=file.path("09_data.zip"))
datadir="09_data"
unzip("09_data.zip",exdir = datadir)
```
## Working with _raw_ HDF files
Will only work if your `gdal` was compiled with HDF support
```{r, eval=T}
gdalinfo(formats = T) %>% grep(pattern="HDF",value=T)
```
```{r}
hdf=file.path(datadir,"MCD12Q1.A2012001.h12v04.051.2014288200441_subset.hdf")
```
Use `gdalinfo()` to print information about the file.
```{r, eval=F}
gdalinfo(hdf)
```
### Subdatasets
An important component of the metadata of a HDF file is the list of 'subdatasets' that are inside the file. HDF files can hold any number of different datasets and you need to use the specific subdataset
* SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:\"09_data/MCD12Q1.A2012001.h12v04.051.2014288200441_subset.hdf\":MOD12Q1:Land_Cover_Type_1
* SUBDATASET_1_DESC=[2400x2400] Land_Cover_Type_1 MOD12Q1 (8-bit unsigned integer)
* SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:\"09_data/MCD12Q1.A2012001.h12v04.051.2014288200441_subset.hdf\":MOD12Q1:Land_Cover_Type_1_Assessment
* SUBDATASET_2_DESC=[2400x2400] Land_Cover_Type_1_Assessment MOD12Q1 (8-bit unsigned integer)
* SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:\"09_data/MCD12Q1.A2012001.h12v04.051.2014288200441_subset.hdf\":MOD12Q1:Land_Cover_Type_QC.Num_QC_Words_01
* SUBDATASET_3_DESC=[2400x2400] Land_Cover_Type_QC.Num_QC_Words_01 MOD12Q1 (8-bit unsigned integer)
#### Translate to GEOtif
```{r, eval=F}
gdal_translate("HDF4_EOS:EOS_GRID:\"09_data/MCD12Q1.A2012001.h12v04.051.2014288200441_subset.hdf\":MOD12Q1:Land_Cover_Type_1",
"test.tif")
gdalinfo("test.tif",nomd=T)
```
#### Plot it
```{r}
d=raster("test.tif")
plot(d)
```
See also the `ModisDownload()` function in `library(rts)`:
* Downloads series of MODIS images in a specific timeframe for specified tile(s)
* MODIS Reproject Tool (MRT) software to mosaic, reproject, reformat
# Use MODISTools package to access the MODISweb
## List MODIS products
```{r}
GetProducts()
```
```{r}
GetBands(Product = "MCD12Q1")
```
## Selection locations
```{r}
loc=rbind.data.frame(
list("UB Spine",43.000753, -78.788195))
colnames(loc)=c("loc","lat","long")
coordinates(loc)=cbind(loc$long,loc$lat)
```
## Available dates
```{r}
mdates=GetDates(Product = "MOD11A2", Lat = loc$lat[1], Long = loc$long[1])
```
### MODIS date codes:
`.A2006001` - Julian Date of Acquisition (A-YYYYDDD)
Convert to a _proper_ date:
* Drop the "`A`"
* Specify date format with julian day `[1,365]`
```{r}
td=mdates[1:5]
td
```
`sub()` to _substitute_ a character in a `vector()`
```{r}
sub("A","",td)
```
Check `?strptime` for date formats.
* `%Y` 4-digit year
* `%j` 3-digit Julian day
```{r}
sub("A","",td)%>%
as.Date("%Y%j")
```
## Add start and end dates to `loc` object
```{r}
dates=mdates%>%sub(pattern="A",replacement="")%>%as.Date("%Y%j")
loc$start.date <- min(as.numeric(format(dates,"%Y")))
loc$end.date <- max(as.numeric(format(dates,"%Y")))
```
## Identify (and create) download folders
Today we'll work with:
* Land Surface Temperature (`lst`): MOD11A2
* Land Cover (`lc`): MCD12Q1
```{r}
lstdir=file.path(datadir,"lst")
if(!file.exists(lstdir)) dir.create(lstdir)
lcdir=file.path(datadir,"lc")
if(!file.exists(lcdir)) dir.create(lcdir)
```
## Download subset
`Size` whole km (integers) for each direction.
`Size=c(1,1)` for 250m resolution data will return a 9x9 pixel tile for each location, centred on the input coordinate.
`Size=c(0,0)` only the central pixel.
**Maximum** size tile `Size=c(100,100)`
This can take a few minutes to run, so you can use the file provided in the data folder.
### Get Land Surface Temperature Data
```{r, eval=F}
MODISSubsets(LoadDat = loc,
Products = c("MOD11A2"),
Bands = c( "LST_Day_1km", "QC_Day"),
Size = c(10,10),
SaveDir=lstdir,
StartDate=T)
```
### Get LULC
```{r, eval=F}
MODISSubsets(LoadDat = loc,
Products = c("MCD12Q1"),
Bands = c( "Land_Cover_Type_1"),
Size = c(10,10),
SaveDir=lcdir,
StartDate=T)
```
List available files:
```{r}
lst_files=list.files(lstdir,pattern="Lat.*asc",full=T)
head(lst_files)
```
Output:
* 1 file per location in `loc`
* Rows: time-steps
* Columns: data bands
```{r}
lst_subset <- read.csv(lst_files[1],header = FALSE, as.is = TRUE)
dim(lst_subset)
lst_subset[1:5,1:15]
```
## Convert to ASCII Grid raster files
Use `MODISGrid()` to convert to separate [ASCII Grid format](http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/spatial_analyst_tools/esri_ascii_raster_format.htm) files:
```
NCOLS xxx
NROWS xxx
XLLCENTER xxx | XLLCORNER xxx
YLLCENTER xxx | YLLCORNER xxx
CELLSIZE xxx
NODATA_VALUE xxx
row 1
row 2
...
row n
```
## Convert LST Data
```{r, eval=F}
MODISGrid(Dir = lstdir,
DirName = "modgrid",
SubDir = TRUE,
NoDataValues=
list("MOD11A2" = c("LST_Day_1km" = 0,
"QC_Day" = -1)))
```
## Convert LandCover Data
```{r, eval=F}
MODISGrid(Dir = lcdir,
DirName = "modgrid",
SubDir = TRUE,
NoDataValues=
list("MCD12Q1" = c("Land_Cover_Type_1" = 255)))
```
## Get lists of `.asc` files
```{r}
lst_files=list.files(file.path(lstdir,"modgrid"),recursive=T,
pattern="LST_Day.*asc",full=T)
head(lst_files)
lstqc_files=list.files(file.path(lstdir,"modgrid"),recursive=T,
pattern="QC_Day.*asc",full=T)
```
## Create raster stacks of evi and evi qc data
```{r}
lst=stack(lst_files)
plot(lst[[1:2]])
```
### Check gain and offset in [metadata](https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod11a2).
```{r}
gain(lst)=0.02
offs(lst)=-273.15
plot(lst[[1:2]])
```
# MODLAND Quality control
See a detailed explaination [here](https://lpdaac.usgs.gov/sites/default/files/public/modis/docs/MODIS_LP_QA_Tutorial-1b.pdf). Some code below from [Steven Mosher's blog](https://stevemosher.wordpress.com/2012/12/05/modis-qc-bits/).
## MOD11A2 (Land Surface Temperature) Quality Control
[MOD11A2 QC Layer table](https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod11a2)
```{r}
lstqc=stack(lstqc_files)
plot(lstqc[[1:2]])
```
### LST QC data
QC data are encoded in 8-bit 'words' to compress information.
```{r}
values(lstqc[[1:2]])%>%table()
```
```{r}
intToBits(65)
intToBits(65)[1:8]
as.integer(intToBits(65)[1:8])
```
#### MODIS QC data are _Big Endian_
Format Digits value sum
---- ---- ---- ----
Little Endian 1 0 0 0 0 0 1 0 65 2^0 + 2^6
Big Endian 0 1 0 0 0 0 0 1 65 2^6 + 2^0
Reverse the digits with `rev()` and compare with QC table above.
```{r}
rev(as.integer(intToBits(65)[1:8]))
```
QC for value `65`:
* LST produced, other quality, recommend exampination of more detailed QA
* good data quality of L1B in 7 TIR bands
* average emissivity error <= 0.01
* Average LST error <= 2K
## Your turn
What does a QC value of 81 represent?
```{r, purl=F}
rev(as.integer(intToBits(81)[1:8]))
# LST produced, other quality, recommend exampination of more detailed QA
# Other quality data
# Average emissivity error <= 0.01
# Average LST error <= 2K
```
### Filter the the lst data using the QC data
```{r}
## set up data frame to hold all combinations
QC_Data <- data.frame(Integer_Value = 0:255,
Bit7 = NA, Bit6 = NA, Bit5 = NA, Bit4 = NA,
Bit3 = NA, Bit2 = NA, Bit1 = NA, Bit0 = NA,
QA_word1 = NA, QA_word2 = NA, QA_word3 = NA,
QA_word4 = NA)
##
for(i in QC_Data$Integer_Value){
AsInt <- as.integer(intToBits(i)[1:8])
QC_Data[i+1,2:9]<- AsInt[8:1]
}
QC_Data$QA_word1[QC_Data$Bit1 == 0 & QC_Data$Bit0==0] <- "LST GOOD"
QC_Data$QA_word1[QC_Data$Bit1 == 0 & QC_Data$Bit0==1] <- "LST Produced,Other Quality"
QC_Data$QA_word1[QC_Data$Bit1 == 1 & QC_Data$Bit0==0] <- "No Pixel,clouds"
QC_Data$QA_word1[QC_Data$Bit1 == 1 & QC_Data$Bit0==1] <- "No Pixel, Other QA"
QC_Data$QA_word2[QC_Data$Bit3 == 0 & QC_Data$Bit2==0] <- "Good Data"
QC_Data$QA_word2[QC_Data$Bit3 == 0 & QC_Data$Bit2==1] <- "Other Quality"
QC_Data$QA_word2[QC_Data$Bit3 == 1 & QC_Data$Bit2==0] <- "TBD"
QC_Data$QA_word2[QC_Data$Bit3 == 1 & QC_Data$Bit2==1] <- "TBD"
QC_Data$QA_word3[QC_Data$Bit5 == 0 & QC_Data$Bit4==0] <- "Emiss Error <= .01"
QC_Data$QA_word3[QC_Data$Bit5 == 0 & QC_Data$Bit4==1] <- "Emiss Err >.01 <=.02"
QC_Data$QA_word3[QC_Data$Bit5 == 1 & QC_Data$Bit4==0] <- "Emiss Err >.02 <=.04"
QC_Data$QA_word3[QC_Data$Bit5 == 1 & QC_Data$Bit4==1] <- "Emiss Err > .04"
QC_Data$QA_word4[QC_Data$Bit7 == 0 & QC_Data$Bit6==0] <- "LST Err <= 1"
QC_Data$QA_word4[QC_Data$Bit7 == 0 & QC_Data$Bit6==1] <- "LST Err > 2 LST Err <= 3"
QC_Data$QA_word4[QC_Data$Bit7 == 1 & QC_Data$Bit6==0] <- "LST Err > 1 LST Err <= 2"
QC_Data$QA_word4[QC_Data$Bit7 == 1 & QC_Data$Bit6==1] <- "LST Err > 4"
kable(head(QC_Data))
```
### Select which QC Levels to keep
```{r}
keep=QC_Data[QC_Data$Bit1 == 0,]
keepvals=unique(keep$Integer_Value)
keepvals
```
### How many observations will be dropped?
```{r,warning=F}
qcvals=table(values(lstqc)) # this takes a minute or two
QC_Data%>%
dplyr::select(everything(),-contains("Bit"))%>%
mutate(Var1=as.character(Integer_Value),
keep=Integer_Value%in%keepvals)%>%
inner_join(data.frame(qcvals))%>%
kable()
```
Do you want to update the values you are keeping?
### Filter the LST Data keeping only `keepvals`
These steps take a couple minutes.
Make logical flag to use for mask
```{r}
lstkeep=calc(lstqc,function(x) x%in%keepvals)
```
Plot the mask
```{r,fig.height=12}
gplot(lstkeep[[4:8]])+
geom_raster(aes(fill=as.factor(value)))+
facet_grid(variable~.)+
scale_fill_manual(values=c("blue","red"),name="Keep")+
coord_equal()+
theme(legend.position = "bottom")
```
Mask the lst data using the QC data
```{r}
lst2=mask(lst,mask=lstkeep,maskval=0)
```
## Add Dates to Z dimension
```{r}
tdates=names(lst)%>%
sub(pattern=".*_A",replacement="")%>%
as.Date("%Y%j")
names(lst2)=1:nlayers(lst2)
lst2=setZ(lst2,tdates)
```
## Summarize to Seasonal climatologies
Use `stackApply()` with a seasonal index.
```{r}
tseas=as.numeric(sub("Q","",quarters(getZ(lst2))))
tseas[1:20]
lst_seas=stackApply(lst2,
indices = tseas,
mean,na.rm=T)
names(lst_seas)=c("Q1_Winter",
"Q2_Spring",
"Q3_Summer",
"Q4_Fall")
```
```{r,fig.height=9}
gplot(lst_seas)+geom_raster(aes(fill=value))+
facet_wrap(~variable)+
scale_fill_gradientn(colours=c("blue",mid="grey","red"))+
coord_equal()+
theme(axis.text.x=element_text(angle=60, hjust=1))
```
## Your turn
Use `stackApply()` to generate and plot monthly median lst values.