BRAG Meeting – 16 May 2013

Geospatial Visualisation (Su & Matt)

Su Yun Kang and Matt Moores presented a seminar in our series on visualisation. The topic was visualisation of geospatial data in R.

Shapefiles

Some simple maps:

library(maps)
library(mapdata)
map("worldHires", "Canada", xlim = c(-141, -53), ylim = c(40, 85),
     col = "brown", fill = TRUE)

canada

map("world", "Australia", col = "blue", fill = TRUE)

australia

map("world","Australia", xlim=c(137.9960, 153.552171),
     ylim=c(-29.1779, -9.142176), col="green", fill=TRUE)

qld
You can also read in your own shapefile:

library(maptools)
library(rworldmap)

region <- readShapeSpatial("SLA_QLD_06_region_new.shp")

mapParams <- mapPolys(mapToPlot=region, nameColumnToPlot="SLA_CODE",
 catMethod="categorical",colourPalette="topo", addLegend="FALSE",
 borderCol='black', mapTitle='Queensland', missingCountryCol=NA,
 xlim=c(137.9960, 153.552171), ylim=c(-29.1779, -9.142176))

QLD

mapParams <- mapPolys(mapToPlot=region, nameColumnToPlot="CATEGORY",
 numCats=6, catMethod="categorical", colourPalette="rainbow",
 addLegend="FALSE", borderCol='black', mapTitle='Queensland',
 missingCountryCol=NA, xlim=c(137.9960, 153.552171),
 ylim=c(-29.1779, -9.142176))

do.call(addMapLegendBoxes, c(mapParams,
 list(legendText=c("Region 1", "Region 2", "Region 3", "Region 4",
 "Region 5"), x='topright', title='', horiz=FALSE)))

QLD2

googleVis

visualise data in a web browser using Google Visualisation API.
Example showing US data by state:

library(googleVis)
require(datasets)
states <- data.frame(state.name, state.x77)
G3 <- gvisGeoMap(states, "state.name", "Illiteracy", 
    options = list(region = "US", dataMode = "regions",
    width = 600, height = 400))
plot(G3)

Example: plot country-level data

data(Exports)
View(Exports)
Geo <- gvisGeoMap(Exports, locationvar = "Country", numvar = "Profit",
    options = list(height = 400, dataMode = "regions"))
plot(Geo)

World population

WorldPopulation = data.frame(Country = Population$Country,
    Population.in.millions = round(Population$Population/1e+06, 0),
    Rank = paste(Population$Country, "Rank:", Population$Rank))
G5 <- gvisGeoMap(WorldPopulation, "Country", "Population.in.millions",
   "Rank", 
    options = list(dataMode = "regions", width = 600, height = 300))
plot(G5)

Example: Plotting point data onto a google map (internet)

data(Andrew)
M1 <- gvisMap(Andrew, "LatLong", "Tip", options = list(showTip = TRUE,
    showLine = F, enableScrollWheel = TRUE, mapType = "satellite",
    useMapTypeControl = TRUE, width = 800, height = 400))
plot(M1)

RgoogleMaps

get maps from Google

library(RgoogleMaps)
newmap <- GetMap(center=c(36.7,-5.9), zoom=8, destfile="newmap.png",
                 maptype = "satellite")
img <- readPNG('newmap.png',native=TRUE)

latRng <- c(newmap$BBOX$ll[1],newmap$BBOX$ur[1])
lonRng <- c(newmap$BBOX$ll[2],newmap$BBOX$ur[2])
plot(lonRng, latRng, type='n', xlab="longitude (degrees)",
     ylab="latitude (degrees)", main="Gibraltar")
rasterImage(img, lonRng[1], latRng[1], lonRng[2], latRng[2])

Gibraltar

try different map types:

newmap2 <- GetMap.bbox(lonR = c(-5, -6), latR = c(36, 37),
    destfile = "newmap2.png", maptype = "terrain")
img <- readPNG("newmap2.png", native = TRUE)
plot(c(-5, -6), c(36, 37), type = "n", xlab = "longitude (degrees)",
    ylab = "latitude (degrees)", main = "")
rasterImage(img, -6, 36, -5, 37)

gibraltar-2

and plot data onto these maps, e.g. these 3 points

PlotOnStaticMap(lat = c(36.3, 35.8, 36.4), lon = c(-5.5, -5.6, -5.8),
    zoom = 10, cex = 2, pch = 19, col = "red", FUN = points, add = F)

gibraltar-3

rasterVis

The 2006 National Land Cover Dataset 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
levelplot(patch)
rasterVis levelplot

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

Yosemite

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, at=my.labs.at)))

NLCD Yosemite

Google Earth

generate KML for Google Earth:

require(R2G2)
require(dismo)
require(XML)
require(gtools)

Overlay2GE(ext[4:1],image="NLCD.png",goo="figure/Overlay2GE.kml")

Links

Steve Mosher
Oscar Perpiñán Lamigueiro
Displaying time series, spatial and space-time data with R: stories of space and time
R2G2 on Recology
Google Earth Engine
googleVis in WordPress
MapBox

Advertisements

One thought on “BRAG Meeting – 16 May 2013

  1. Pingback: Posterior samples | Sam Clifford

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