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