### Advanced Spatial Methods in R

```ADVANCED SPATIAL
METHODS IN R
Michael Mann
George Washington University
Department of Geography
[email protected]
http://michaelmann.i234.me/
Overview
• Review Basics
• Setting up Space-time data
• Space-time plots
• library(RasterVis)
• library(plotKML)
• Vector & Other Data Visualization
• library(ggplot2)
• Mapping ggplot2 visualizations
• library(ggmap)
Before we begin
Best places for help
• http://stackoverflow.com/
• Question & Answer form. Quick & high quality
responses
• ?function_name
• Look up help files for a function from any library
• http://gis.stackexchange.com/
• Similar to stackoverflow, but targeted to spatial
community
• The fridge – grab a beer and spend some time
Interpreting My Slides!
Inputs into command line
s = c("aa", "bb", "cc", "dd", "ee")
s
Outputs
[1] "aa” "bb” "cc” "dd” "ee"
Notes
An important
note
Data is in the MMSpatialData folder
Quick survey
• Please raise you hand if (before
today) you have never used R
or a similar language.
Let’s even the playing field
Beginners
• Please look around you. Move if
there is a beginner cluster!
Experts
• Put on your teaching hats!
Remember how difficult this
material is. Make sure to help
Overview
Indexing a vector
Indexing a data.frame
s = c("aa", "bb", "cc", "dd", "ee")
s = data.frame(col1=c(1,2,3),
col2=c(5,6,7))
s[1]
s[c(2,5)]
[1] ‘aa’
[1] ‘bb’ ‘ee’
col1 col2
1 1 5
2 2 6
3 3 7
s[c(2:5)]
s[2,1]
s[,’col2’]
[1] ‘bb’ ‘cc’ ‘dd’ ‘ee’
[1] 2
[1] 5 6 7
vector_name[ position_# ]
vector_name[ row#, col# ]
Objective 1 – Raster space-time plots
• library(raster)
• Raster Stacks & Bricks
• Multidimensional raster objects
• Multi-layer (red,green,blue)
• Multi-dim (time series, multi variable)
stack[row# ,col#, Z#]
Y
Z
X
Objective 1 – Raster space-time plots
• Multi-layer Raster ‘Brick’
b <- brick(system.file("external/rlogo.grd", package="raster)
plot(b)
Brick[X,Y,Z]
Objective 1 – Raster space-time plots
• Multi-layer Raster Brick
plotRGB(b, r=1, g=2, b=3)
Objective 1 – Raster space-time plots
• Multi-dimentional Raster Stack
• Good for time-series of rasters, or multivariate analysis
Time Series
stack[x,y, c(cwd2000-2010) ]
Multivariate
stack[x,y, c(crime,pop) ]
Objective 1 – Raster space-time plots
• Multi-dimentional Raster Stack
• Good for time-series of rasters, or multivariate analysis
Time Series
stack[x,y, c(cwd2000-2010) ]
Use: Data
Visualization
Multivariate
stack[x,y, c(crime,pop) ]
Use: Regression
Analysis, modeling
Task 1.1 – Setup Raster Stack Data
grep()
dir()
dir(‘C://SESYNC//data’)
grep(‘a’,
[1]
grep(‘c’,
[1]
[1] ‘data.zip’ ‘ggmap vinette.pdf’….
c( ‘a’, ‘b’, ‘c’, ‘a’ )
1
4
c( ‘tab’, ‘car’, ‘bat’ )
2
)
)
Task 1.1 – Setup Raster Stack Data
Helpful Functions – Using grep to index a vector
Find the location of a element with ‘c’ in it
s = c("aa", "bb", "cc", "dd", "ee")
grep(‘c’, s )
[1]
3
Query vector s with the grep
s[ grep(‘c’, s ) ]
[1]
‘cc’
Task 1.1 – Setup Raster Stack Data
Helpful Functions – Using grep to index a vector
Find the location of a element with ‘c’ in it
s = c("aa", "bb", "cc", "dd", "ee")
grep(‘c’, s )
[1]
3
Query vector s with the grep
s[ grep(‘c’, s ) ]
[1]
‘cc’
Task 1.1 – Let’s get started!
Task 1.2 – Create Raster Stacks and
Assign Name Labels
paste()
paste(‘Hi',c('Bill','Bob', 'Sam'), sep=' ')
names()
names(rstack) = c(‘test2001’,’test2002’)
[1] “Hi Bill" “Hi Bob" “Hi Sam"
paste(’aet',c(’2001',’2002’), sep=’_')
[1] "aet_2001" "aet_2002"
test2002
test2001
Task 1.2 – Create Raster Stacks and
Assign Name Labels
paste()
paste(‘Hi',c('Bill','Bob', 'Sam'), sep=' ')
names()
names(rstack) = c(‘test2001’,’test2002’)
[1] “Hi Bill" “Hi Bob" “Hi Sam"
paste(’aet',c(’2001',’2002’), sep=’_')
[1] "aet_2001" "aet_2002"
test2002
test2001
Task 1.2 – Create Raster Stacks and
Assign Time Labels
seq() & as.Date()
seq(1,15, by=3)
[1] 1 4 7 10 13
Years = seq(as.Date(’2001-01-01'),
as.Date('2010-01-01'), by=’1 year')
[1] "2000-01-01" "2001-01-01"
…
[10] "2009-01-01" "2010-01-01"
setz()
raster_stack = setZ(raster_stack, years)
Task 1.2 – Create Raster Stacks and
Assign Time Labels
seq() & as.Date()
seq(1,15, by=3)
[1] 1 4 7 10 13
Years = seq(as.Date(’2001-01-01'),
as.Date('2010-01-01'), by=’1 year')
[1] "2000-01-01" "2001-01-01"
…
[10] "2009-01-01" "2010-01-01"
setz()
raster_stack = setZ(raster_stack, years)
Important: setZ must be
passed a series of ‘Dates’
(created with as.Date
function)
Task 1.2 – Let’s get started!
Task 1.3 – Visualize Stack Data
Method 1
plot(raster_stack[[1]]) # plot first raster
plot(raster_stack[[2]]) # plot second raster
NOTE: [[ ]] is used b/c stack is a list object
Method 2
plot( raster(raster_stack,'HDen_1989') )
NOTE: raster() is used b/c… well that is just
how it works
Task 1.3 – Let’s get started!
rasterVis & plotKML packages
rasterVis
• Excellent tutorial available at http://rastervis.r-forge.r-project.org/
• Data Format
• Raster stack with z-dim set to dates using: as.Date()
• Spatial points or polygons data.frame with z-dim set to dates
Hovmoller Plot
Horizon Plot
rasterVis & plotKML packages
rasterVis
• Data Format 2
• Raster stack OR Spatial points or polygons data.frame with zdim set to slope, or direction
Vectorplot
Stream Plot
rasterVis & plotKML packages
plotKML
• Excellent tutorial: http://plotkml.r-forge.r-project.org/
• Data Format
• Exports many formats including raster stacks
Note: This code outputs both a
Housing.kml file and a series of other
image files (.png files). In order for this
to work, the Housing.kml file needs to
be in the same directory as all the
image files.
rasterVis & plotKML packages
plotKML
• Excellent tutorial: http://plotkml.r-forge.r-project.org/
• Data Format
• Exports many formats including raster stacks
All data must be unprojected Lat Lon
"+proj=longlat +datum=WGS84”
Note: This code outputs both a
Housing.kml file and a series of other
image files (.png files). In order for this
to work, the Housing.kml file needs to
be in the same directory as all the
image files.
Task 1.5 – Let’s get started!
Task 1.6 – Let’s get started!
Objective 2: Intro to ggplot2
• One of the best data visualization tools in R
• Documentation available here: http://ggplot2.org/
Objective 2: Intro to ggplot2
• A plot is made up of multiple layers.
• A layer consists of data, a set of mappings between
variables and aesthetics, a geometric object and a
statistical transformation
• Scales control the details of the mapping.
• All components are independent and reusable.
Objective 2: Intro to ggplot2
• plot
• aesthetics
• geometric object
• scales control
Typical Command
a = ggplot(movies, aes(y = budget, x = year, group = round_any(year, 10) ))
+ geom_boxplot() + scale_y_log10()
plot(a)
ggplot2 uses data.frames!
aggregate()
Summarizes data of interest by factors (categorical data)
aggregate( rating ~ year ,data= movies,
FUN='mean')
fortify()
Converts spatial polygons, lines, points to data.frame
usable in ggplot2
jepson.points = fortify(jepson,
region="id")
Objective 2: Intro to ggplot2
• plot reminder
• aesthetics
• geometric object
• scales control
Typical Command
a = ggplot(jepson.df) + aes(long,lat,group=group,fill=ECOREGION) + geom_polygon()
plot(a)
Objective 3: Intro to ggmap
• Ggmap enables visualization of layered graphics using
implementation similar to ggplot2
• Combines the functionality of ggplot2 and spatial
information of static maps from Google Maps,
Objective 3: Intro to ggmap
• Ggmap enables visualization of layered graphics using
implementation similar to ggplot2
• Defined by a central Lat and Lon and a ‘Zoom’ level
2. Takes additional ‘+’ commands to overlay ggplot2
graphics
Zoom = 5
Zoom = 14
Objective 3: Intro to ggmap
• plot
• aesthetics
• geometric object
• scales control
• coordinate system
• Done in background through qmap
Typical Command
HoustonMap <- qmap("houston", zoom = 13, color = "bw")
HoustonMap + geom_point(data=violent_crimes,aes(x = lon, y = lat, colour = offense ) )
plot( HoustonMap )
projectRaster() & project()
• Converts spatial -polygons, -lines, -points to data.frame
• usable in ggplot2 & ggmap.
• Ggmap data must be in unprojected lat lon (defined below)
# Project a raster to unprojected Lat Lon using nearest neighbor algorithm
stack_proj = projectRaster(raster_stack, crs="+proj=longlat +datum=WGS84”
, method='ngb')
# Project a polygon to unprojected Lat Lon using nearest neighbor algorithm
jepson_proj = project(jepson, crs="+proj=longlat +datum=WGS84”)
fortify()
Converts spatial polygons, lines, points to data.frame
usable in ggplot2
jepson.points = fortify(jepson, region="id")
geocode()
Converts text
names to Lat Lon
coordindates
geocode("the white house")
revgeocode()
Converts Lat Lon to text
revgeocode(c(-77.03650, 38.89768),
subset()
data = data.frame(name=c(‘mike’,’john’,’jim’), age=c(4,3,6))
1
2
3
name age
mike 4
john
3
jim
6
subset(data, name != ‘mike’ )
subset(data, age > 3 )
2
3
name age
john
3
jim
6
1
3
name age
mike 4
jim
6
Task 3.1: Learn Basic ggmap Functions
Task 3.2: Crime Mapping In Houston TX