title: "Satellite Remote Sensing"

```{r, echo=FALSE, message=FALSE, results='hide', purl=FALSE}

<div>

<iframe src="09_presentation/09_RemoteSensing.html" width="75%" height="400px"> </iframe>

</div>

[<i class="fa fa-file-code-o fa-3x" aria-hidden="true"></i> 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.

Will only work if your `gdal` was compiled with HDF support

```{r, eval=T}

Use `gdalinfo()` to print information about the file.

```{r, eval=F}

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)

```{r, eval=F}

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

`.A2006001` - Julian Date of Acquisition (A-YYYYDDD)

Convert to a _proper_ date:

* Drop the "`A`"

* Specify date format with julian day `[1,365]`

`sub()` to _substitute_ a character in a `vector()`

Check `?strptime` for date formats.

* `%Y` 4-digit year

* `%j` 3-digit Julian day

Today we'll work with:

* Land Surface Temperature (`lst`): MOD11A2

* Land Cover (`lc`): MCD12Q1

`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.

```{r, eval=F}

```{r, eval=F}

List available files:

Output:

* 1 file per location in `loc`

* Rows: time-steps

* Columns: data bands

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:

```{r, eval=F}

```{r, eval=F}

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 QC Layer table](https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod11a2)

QC data are encoded in 8-bit 'words' to compress information.

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.

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

<div class="well">

What does a QC value of 81 represent?

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo1">Show Solution</button>

<div id="demo1" class="collapse">

```{r, purl=F}

</div>

</div>

Do you want to update the values you are keeping?

These steps take a couple minutes.

Make logical flag to use for mask

Plot the mask

Mask the lst data using the QC data

Use `stackApply()` with a seasonal index.

<div class="well">

Use `stackApply()` to generate and plot monthly median lst values.

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo2">Show Solution</button>

<div id="demo2" class="collapse">

```{r, purl=F}

</div>

</div>

```{r, warning=F}

See the `library(rts)` for more timeseries related functions.

Get cover clases from [MODIS website](https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd12q1)

```

Convert to `factor` raster

Warnings about `.checkLevels()` OK here because some factors not present in this subset...

Extract values from `lst` and `lc` rasters.

Melt table and add LandCover Name

<div class="well">

Calculate the maximum observed seasonal average lst in each land cover type.

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo3">Show Solution</button>

<div id="demo3" class="collapse">

```{r, purl=F}

</div>

</div>

Things to think about:

* What tests would you use to identify differences?

* Do you need to worry about unequal sample sizes?