---
title: "Landsat 8 Vegetation Index"
author: "Matt Moores"
date: "15/03/2017"
output: html_document
---

Region of interest is the city of Brisbane, Australia:
```{r}
library(sp)
library(rgdal)
library(ggmap)

latlon <- matrix(c(-27.466667, 153.033333), nrow=1)
colnames(latlon) <- c("lat","lon")
coords <- data.frame(latlon)
coordinates(coords) <- c("lon","lat")
proj4string(coords) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
img <- get_googlemap(c(latlon[1,2], latlon[1,1]), zoom=9, maptype="terrain", filename="data/Brisbane_GGmap")
plot(img)
```

Read the GeoTIFF file into a raster object and crop to 30km square:
```{r}
library(raster)
library(rasterVis)

fileName <- "LC80890792016334LGN00" # Nov 29, 2016
#ndvi_rstr <- raster(paste0("data/",fileName,"_sr_ndvi.tif"))
#extent(ndvi_rstr)
#proj4string(ndvi_rstr)
#coords <-  spTransform(coords, proj4string(ndvi_rstr))
#Lon <- coordinates(coords)[1,"lon"]
#Lat <- coordinates(coords)[1,"lat"]
#e <- extent(Lon-15000,Lon+15000,Lat-17000,Lat+13000) # 30km square ROI
#system.time(patch <- crop(ndvi_rstr,e))
#levelplot(patch/10000, par.settings=viridisTheme())
```

Read metadata and plot histogram:
```{r}
library(RStoolbox)
#usgs_mtl <- readMeta(paste0("data/",fileName,"_MTL.txt"))
#tm_mtl <- as.Date(usgs_mtl$ACQUISITION_DATE)
#newName <- format(tm_mtl, "%Y%m%d")

#ndvi <- as.matrix(patch)/10000
load("data/NDVIp089r079_20161129.rda")
dim(ndvi)
sum(is.na(ndvi))
range(ndvi)
hist(ndvi, xlim=c(-1,1), col=3, main="", xlab="NDVI", freq=FALSE, breaks=50)
lines(density(ndvi), col=4, lty=2, lwd=2)
```

If plotting a matrix, need to transpose and flip the Y axis:
```{r}
levelplot(t(ndvi)[,nrow(ndvi):1], par.settings=viridisTheme())
```

