---
title: "MODIS Ocean Color"
author: "Matt Moores"
date: "15/03/2017"
output: html_document
---

Get MODIS tile ID for a given latitude and longitude:
```{r}
library(sp)
library(raster)
library(maps)
library(mapdata)
library(MODIS)

y0 <- -10
x0 <- 142.3
y1 <- -20
x1 <- 149
tileID <- getTile(extent=list(xmin=x0, xmax=x1, ymin=y0, ymax=y1))
print(tileID$tile)
map("worldHires","Australia", xlim=c(x0, x1), ylim=c(y1, y0), col="red")
map.cities(country="Australia", label=TRUE, minpop=10000, cex=1, lwd=2) # Cairns & Townsville
```

Read the GeoTIFF file into a raster object and plot:
```{r}
chl <- raster("data/A2013320041000_reprojected.tif")
chl[chl==-32767] <- NA
plot(chl)
projection(chl)
```

Crop to the region of interest, far north Queensland:
```{r}
library(rasterVis)

mycorners <- SpatialPoints(list(x=c(x0,x1),y=c(y0,y1)))
proj4string(mycorners) <- "+init=epsg:4326"
mycornersm6 <- spTransform(mycorners, CRS = CRS(projection(chl)))
chl_sub <- crop(chl, extent(mycornersm6))
levelplot(chl_sub)
```

Plot the coastline and label some cities:
```{r}
plot(chl_sub, col = topo.colors(341), xlab="longitude", ylab="latitude")
map("worldHires","Australia", xlim=c(x0, x1), ylim=c(y1, y0), col="red", add=TRUE)
map.cities(country="Australia", label=TRUE, minpop=10000, cex=1.5, lwd=2) # Cairns & Townsville
```

Histogram of chlorophyll concentration (with kernel density estimate):
```{r}
hist(chl_sub,breaks=50,freq=FALSE,col="lightgreen",xlab=expression(paste("Chlorophyll a ", (mg/m^3))), main="")
lines(density(as.vector(chl_sub), na.rm=TRUE), lwd=3, col=4)
```
