---
title: "Land Cover"
author: "Matt Moores"
date: "15/03/2017"
output: html_document
---

Read satellite image as a raster object:
```{r}
library(sp)
library(raster)
nlcd <- raster("~/Downloads/nlcd_2011_landcover_2011_edition_2014_10_10/nlcd_2011_landcover_2011_edition_2014_10_10.img")
nlcd_proj <- projection(nlcd)
nlcd_ext <- extent(nlcd)

library(maps)
library(mapdata)
ext_df <- data.frame(x=c(nlcd_ext@xmin, nlcd_ext@xmax), y=c(nlcd_ext@ymin, nlcd_ext@ymax))
coordinates(ext_df) <- c("x", "y") # class SpatialPoints
proj4string(ext_df) <- nlcd_proj
coords <-  spTransform(ext_df, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(coords@bbox[1,], coords@bbox[2,], type='n', xlab="longitude (degrees)",ylab="latitude (degrees)")
map("worldHires", xlim=coords@bbox[1,], ylim=coords@bbox[2,], col=2, add=TRUE)
map.cities(country="USA", label=TRUE, minpop=1e6, cex=0.7, lwd=2)
```

Crop to the region of interest:
```{r}
library(rasterVis)
place <- c("Mariposa", "Half Moon Bay", "San Jose", "U Minnesota")
latlon <- matrix(c(c(37.485, -119.966389),c(37.458889, -122.436944),c(37.333333, -121.9),c(44.975278, -93.234167)),nrow=4,ncol=2,byrow=T)
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")
coords <-  spTransform(coords, CRS(projection(nlcd)))
Lon <- coordinates(coords)[1,"lon"]
Lat <- coordinates(coords)[1,"lat"]
e <- extent(Lon-110000,Lon+110000,Lat-110000,Lat+110000)
system.time(patch <- crop(nlcd,e)) # this takes a few seconds to run
levelplot(patch)
```

```{r}
library(ggmap)
ext <- t(as.matrix(e))
ext[1,] <- ext[1,] - 50000
ext[2,] <- ext[2,] + 50000
ext <- SpatialPoints(ext, CRS(projection(nlcd)))
ext <- coordinates(spTransform(ext, 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/Mariposa_STmap")
plot(ext[,1],ext[,2], type='n',xlab="longitude (degrees)",ylab="latitude (degrees)",main="Yosemite")
rasterImage(img, ext[1,1], ext[1,2], ext[2,1], ext[2,2])
```

```{r}
##   land class codes
LCcodes  <- c(0,11,12,21,22,23,24,31,41,42,43,
              51,52,71,72,73,74,81,82,90,95)
##     a data frame for substituting the landclass values to an  0..N
dex     <- data.frame(id=LCcodes,v=seq(0,20))
  ##    land cover names
  LandCovers <-c(
    "Water",
    "IceSnow",
    "Isa20",
    "Isa50",
    "Isa75",
    "Isa100",
    "Barren",
    "DeciduousForest",
    "EvergreenForest",
    "MixedForest",
    "DwarfScrub",
    "ShrubScrub",
    "GrassHerbaceous",
    "SedgeHerbaceous",
    "Lichens",
    "Moss",
    "PastureHay",
    "CultCrops",
    "WoodyWetlands",
    "EmergentHerbWet" )
lcCol <- c("blue4","white","deeppink","red","red4","grey","grey50","tan",
           "green","green4","springgreen","olivedrab","seagreen","lawngreen",
           "yellow","brown","chocolate","burlywood","yellowgreen","skyblue","blue")
system.time(patch <- subs(patch,dex))
rng <- c(1,20)
## define the breaks of the color key
my.at <- seq(rng[1]-1, rng[2])
## the labels will be placed vertically centered
my.labs.at <- seq(rng[1], rng[2])-0.5
levelplot(patch, at=my.at, margin=FALSE, col.regions=lcCol,
          colorkey=list(labels=list(
            labels=LandCovers, ## classes names as labels
            at=my.labs.at)),sub="National Land Cover Database 2011",main="Land use near Yosemite")
```

```{r}
lcSimpleLabels <- c("water","snow","developed","barren","forest","scrub",
                    "grassland","cultivated","wetlands")
lcSimpleColour <- c("blue","white","red","grey","green4","greenyellow",
                    "green","chocolate1","aquamarine")
LCcodes  <- c(0,11,12,21,22,23,24,31,41,42,43,
              51,52,71,72,73,74,81,82,90,95)
lcSimpleCodes  <- c(0,1,2,3,3,3,3,4,5,5,5,
                    6,6,7,7,7,7,8,8,9,9)
dex <- data.frame(id=LCcodes,v=lcSimpleCodes)
system.time(patch <- subs(patch,dex)) # less than 10min
rng <- c(1,9)
my.at <- seq(rng[1]-1, rng[2])
## the labels will be placed vertically centered
my.labs.at <- seq(rng[1], rng[2])-0.5
#png(filename="figures/simpleNLCD.png",width=1000,height=1000,units="px",bg = "transparent",type = "cairo-png")
levelplot(patch, at=my.at, margin=FALSE, col.regions=lcSimpleColour,
          colorkey=list(labels=list(
            labels=lcSimpleLabels, ## classes names as labels
            at=my.labs.at)))
```

