Advanced Spatial Methods in R

1. Advanced Spatial Methods in R Michael Mann George Washington University Department of Geography mmann1123@gwu.edu http://michaelmann.i234.me/

2. 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)

3. 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

4. Interpreting My Slides! Inputs into command line OutputsNotes Data is in the MMSpatialData folder s = c("aa", "bb", "cc", "dd", "ee") s [1] "aa” "bb” "cc” "dd” "ee" An important note

5. Your location Quick survey • Please raise you hand if (before today) you have never used R or a similar language.

6. 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 your teammates!

7. Overview Helpful Functions Indexing a data.frame Indexing a vector s = data.frame(col1=c(1,2,3), col2=c(5,6,7)) s = c("aa", "bb", "cc", "dd", "ee") s[1] s[c(2,5)] col1 col2 1 1 5 2 2 6 3 3 7 [1] ‘aa’ [1] ‘bb’ ‘ee’ s[2,1] s[,’col2’] s[c(2:5)] [1] 2 [1] 5 6 7 [1] ‘bb’ ‘cc’ ‘dd’ ‘ee’ vector_name[ row#, col# ] vector_name[ position_# ]

8. 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) Indexing Your Data stack[row# ,col#, Z#] Y Z X

9. 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]

10. Objective 1 – Raster space-time plots • Multi-layer Raster Brick plotRGB(b, r=1, g=2, b=3)

11. Objective 1 – Raster space-time plots • Multi-dimentional Raster Stack • Good for time-series of rasters, or multivariate analysis Time Series Multivariate stack[x,y, c(cwd2000-2010) ] stack[x,y, c(crime,pop) ]

12. Objective 1 – Raster space-time plots • Multi-dimentional Raster Stack • Good for time-series of rasters, or multivariate analysis Time Series Multivariate stack[x,y, c(cwd2000-2010) ] stack[x,y, c(crime,pop) ] Use: Regression Analysis, modeling Use: Data Visualization

13. Task 1.1 – Setup Raster Stack Data Helpful Functions grep() dir() grep(‘a’, c( ‘a’, ‘b’, ‘c’, ‘a’ ) ) dir(‘C://SESYNC//data’) [1] 1 4 grep(‘c’, c( ‘tab’, ‘car’, ‘bat’ ) ) [1] 2 [1] ‘data.zip’ ‘ggmapvinette.pdf’….

14. 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 Query vector s with the grep s = c("aa", "bb", "cc", "dd", "ee") grep(‘c’, s ) [1] 3 s[ grep(‘c’, s ) ] [1] ‘cc’

15. 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 Query vector s with the grep s = c("aa", "bb", "cc", "dd", "ee") grep(‘c’, s ) [1] 3 TRY THIS! s[ grep(‘c’, s ) ] [1] ‘cc’

16. Task 1.1 – Let’s get started!

17. Task 1.2 – Create Raster Stacks and Assign Name Labels Helpful Functions names() paste() names(rstack) = c(‘test2001’,’test2002’) paste(‘Hi',c('Bill','Bob', 'Sam'), sep=' ') [1] “Hi Bill" “Hi Bob" “Hi Sam" paste(’aet',c(’2001',’2002’), sep=’_') [1] "aet_2001" "aet_2002" test2002 test2001

18. Task 1.2 – Create Raster Stacks and Assign Name Labels Helpful Functions names() paste() names(rstack) = c(‘test2001’,’test2002’) paste(‘Hi',c('Bill','Bob', 'Sam'), sep=' ') [1] “Hi Bill" “Hi Bob" “Hi Sam" TRY PASTE()! paste(’aet',c(’2001',’2002’), sep=’_') [1] "aet_2001" "aet_2002" test2002 test2001

19. Task 1.2 – Create Raster Stacks and Assign Time Labels Helpful Functions setz() seq() & as.Date() raster_stack = setZ(raster_stack, years) 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"

20. Task 1.2 – Create Raster Stacks and Assign Time Labels Helpful Functions setz() seq() & as.Date() raster_stack = setZ(raster_stack, years) seq(1,15, by=3) [1] 1 4 7 10 13 Important: setZ must be passed a series of ‘Dates’ (created with as.Datefunction) 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"

21. Task 1.2 – Let’s get started!

22. Task 1.3 – Visualize Stack Data Indexing Your 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

23. Task 1.3 – Let’s get started!

24. Task 1.4 – Challenge Questions

25. Task 1.5 Visualize Space-Time Data rasterVis & plotKMLpackages 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

26. Task 1.5 Visualize Space-Time Data rasterVis & plotKMLpackages rasterVis • Data Format 2 • Raster stack OR Spatial points or polygons data.frame with z-dim set to slope, or direction Vectorplot Stream Plot

27. Task 1.5 Visualize Space-Time Data rasterVis & plotKMLpackages 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.

28. Task 1.5 Visualize Space-Time Data rasterVis & plotKMLpackages plotKML • Excellent tutorial: http://plotkml.r-forge.r-project.org/ • Data Format • Exports many formats including raster stacks • All data must be unprojectedLat 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.

29. Task 1.5 – Let’s get started!

30. Task 1.6 – Let’s get started!

31. Objective 2: Intro to ggplot2 • One of the best data visualization tools in R • Documentation available here: http://ggplot2.org/

32. Objective 2: Intro to ggplot2 • A plot is made up of multiple layers. • A layer consists of data, a set of mappingsbetween variables and aesthetics, a geometric object and a statistical transformation • Scales control the details of the mapping. • All components are independent and reusable.

33. Objective 2: Intro to ggplot2 • plot • aesthetics • geometric object • scalescontrol • ‘+’ add mappings Typical Command a = ggplot(movies, aes(y = budget, x = year, group = round_any(year, 10) )) + geom_boxplot() + scale_y_log10() plot(a)

34. Task 2.1 Setting up your data ggplot2 uses data.frames!

35. Task 2.1 Setting up your data Helpful Functions aggregate() Summarizes data of interest by factors (categorical data) aggregate( rating ~ year ,data= movies, FUN='mean')

36. Task 2.1 Setting up your data Helpful Functions fortify() Converts spatial polygons, lines, points to data.frame usable in ggplot2 jepson.points = fortify(jepson, region="id")

37. Task 2.1 Let’s get going!

38. Objective 2: Intro to ggplot2 • plot reminder • aesthetics • geometric object • scalescontrol • ‘+’ add mappings Typical Command a = ggplot(jepson.df) + aes(long,lat,group=group,fill=ECOREGION) + geom_polygon() plot(a)

39. Task 2.2 Let’s get going!

40. Task 2.3 Let’s get going!

41. Task 2.4 Challenge questions

42. 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, OpenStreetMap, Stamen Maps or CloudMade

43. Objective 3: Intro to ggmap • Ggmap enables visualization of layered graphics using implementation similar to ggplot2 • qmap() Downloads static maps from google, or OSM • Defined by a central Lat and Lon and a ‘Zoom’ level • Takes additional ‘+’ commands to overlay ggplot2 graphics Zoom = 14 Zoom = 5

44. Objective 3: Intro to ggmap • plot • aesthetics • geometric object • scalescontrol • coordinate system • Done in background through qmap • ‘+’ add mappings Typical Command HoustonMap <- qmap("houston", zoom = 13, color = "bw") HoustonMap+ geom_point(data=violent_crimes,aes(x = lon, y = lat, colour = offense ) ) plot( HoustonMap )

45. Task 3.1 Setting up your data Helpful Functions projectRaster() & project() • Converts spatial -polygons, -lines, -points to data.frame • usable in ggplot2 & ggmap. • Ggmap data must be in unprojectedlatlon (defined below) # Project a raster to unprojectedLat Lon using nearest neighbor algorithm stack_proj= projectRaster(raster_stack, crs="+proj=longlat +datum=WGS84” , method='ngb') # Project a polygon to unprojectedLat Lon using nearest neighbor algorithm jepson_proj= project(jepson, crs="+proj=longlat +datum=WGS84”)

46. Task 3.1 Setting up your data Helpful Functions fortify() Converts spatial polygons, lines, points to data.frame usable in ggplot2 jepson.points = fortify(jepson, region="id")

47. Task 3.1 Setting up your data Helpful Functions geocode()revgeocode() Converts text addresses or location names to Lat Lon coordindates Converts LatLon to text addresses revgeocode(c(-77.03650, 38.89768), output = c("address")) geocode("the whitehouse")

48. Task 3.1 Setting up your data Helpful Functions subset() data = data.frame(name=c(‘mike’,’john’,’jim’), age=c(4,3,6)) name age mike 4 john 3 jim 6 name age 2 john 3 3 jim 6 subset(data, name != ‘mike’ ) name age 1 mike4 3 jim 6 subset(data, age > 3 )

49. Task 3.1: Learn Basic ggmap Functions

50. Task 3.2: Crime Mapping In Houston TX