Skip to content

Land Cover

May 6, 2013

Following up from posts by Steve Mosher and Oscar Perpiñán Lamigueiro, I thought I’d try my hand at representing categorical data in rasterVis.

The 2006 National Land Cover Data is a 16GB image that gives land use in one of several categories for the continental USA, at a spatial resolution of 30m per pixel. It is derived from an unsupervised classification of Landsat ETM+ satellite data.

Rather than looking at the entire dataset, we can crop out a region of interest:

library(raster)
library(rasterVis)
library(rgdal)
library(colorspace)

nlcd <- raster("c:/dev/NLCD/nlcd2006_landcover_4-20-11_se5.img")
projection(nlcd)

[1] "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

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
save(patch, file="data/patch.rda")
levelplot(patch)
heatmap of land use in eastern California

heatmap of land use in eastern California

The region of interest corresponds to a 220km square chunk of eastern California. We can obtain a map of the same region by translating the bounding box back into WGS84 coordinates and calling ggmap:

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_stamenmap(as.vector(t(ext)),filename="data/Mariposa_STmap",zoom=9)
save(img, file="data/Mariposa_STmap.rda")
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])

Stamen map

The heatmap above doesn’t really convey much information, since the NLCD pixel values are categorical, not continuous. What we really want to do is assign a colour to each category, as follows:

##   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 2006",main="Land use near Yosemite")

NLCD Yosemite
Mapping the pixel values to a sequential numbering scheme took around 7 minutes on my 3.33GHz Core i5. You’d definitely want a parallel version of the subs function if you were going to do this a lot. Now the cities of Fresno, Modesto and Stockton can be seen in shades of red, while the Sierra Navadas are in grey with white peaks and Mono Lake is deep blue.

However, the 220km image is still very busy, so I mapped the categories down to a simplified set of 9:

load("data/patch.rda")
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])
my.labs.at <- seq(rng[1], rng[2])-0.5
levelplot(patch, at=my.at, margin=FALSE, col.regions=lcSimpleColour,
          colorkey=list(labels=list(
            labels=lcSimpleLabels, ## classes names as labels
            at=my.labs.at)))

NLCD Yosemite

Advertisements

From → Imaging, R

2 Comments
  1. Steven Mosher permalink

    Very Nice matt!

Trackbacks & Pingbacks

  1. Warwick R User Group | Matt Moores

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Richard Everitt's blog

Computational Bayesian statistics

Let's Look at the Figures

David Firth's blog

Nicholas Tierney

Computational Bayesian statistics

One weiRd tip

Computational Bayesian statistics

Series B'log

discussion blog for JRSS Series B papers

Mad (Data) Scientist

Musings, useful code etc. on R and data science

R-bloggers

R news and tutorials contributed by (750) R bloggers

Another Astrostatistics Blog

The random musings of a reformed astronomer ...

Darren Wilkinson's research blog

Statistics, computing, data science, Bayes, stochastic modelling, systems biology and bioinformatics

CHANCE

Computational Bayesian statistics

StatsLife - Significance magazine

Computational Bayesian statistics

(badness 10000)

Computational Bayesian statistics

Igor Kromin

Computational Bayesian statistics

Statisfaction

I can't get no

Xi'an's Og

an attempt at bloggin, nothing more...

Sam Clifford

Postdoctoral Fellow, Bayesian Statistics, Aerosol Science

Bayesian Research & Applications Group

Frontier Research in Bayesian Methodology & Computation

%d bloggers like this: