Well it’s been long enough since my last post. Had a few things on my plate (vacation, holidays, another holiday, some more holidays, and quite a lot of research). March is almost here but the good news is that I have plenty of work stored up to start serving out some intuitive approaches for learning R. Speaking of that…

In the hefty amounts of research I’ve been doing lately, I’ve come across many, MANY R-based blogs and tutorials. There are so many fantastic resources out there. But there are also a few not-so-good ones. Some code examples seem to confuse more than clarify. Scrolling to the bottom of a tutorial for a glance at the comments is usually a good way to gauge whether or not the audience received it well (I’ve also noticed that R-learners are much less negative and troll-like than the majority of those who comment on say, well, every other community-based website in the world). A couple tutorials don’t even include code. Unfortunately, I think often times these bloggers have ulterior motives (showing off their technical/statistical/whatever capacity), consequently flushing a graspable, empowering learning experience down the toilet.

With this site, I’m going to continue attempting to hammer on what is actually transpiring in R, ideally without dragging my feet and stagnating the more advanced users, using an intuitive approach so readers understand not just THAT something happens, but HOW and WHY it happens. This hopefully means they remember, are able to reproduce results, and ultimately grow in their learning. So please, keep the feedback coming (even if it is in troll form)!

Alright enough about me. Let’s pick up from where we left off.

For this post, I am going to show you how to plot, or overlay, the polygons of a shapefile on top of a Google Map. The polygons in this example will be of neighborhoods in the city of Baltimore.

The City of Baltimore is a kind of paragon when it comes to a municipality’s dissemination of public data. Rob Mealy’s amazing blog (on R and maps n’ stuff) first tipped me off to this website. I downloaded the neighborhood shapefiles (Neighborhood 2010.zip) from the site and unzipped the file to my C drive.

Now since we are going to be reading shapefiles into R, we need to install a package that is capable of doing so. There are several, but for this example we are going to use rgdal.

install.packages('rgdal')
library(rgdal)

Since I unzipped my shapefile data to my C drive, I am going to tell R it is from THERE I will be working. I set this as my working directory with:

setwd("C:/")

From now on during this session, R will automatically use this location to retrieve and save files, unless specifically told to do otherwise.

We read in the shapefile with:

Neighborhoods <- readOGR(".","nhood_2010")

We named the shapefile “Neighborhoods” (by typing the name of to the left of <-). The first set of quotations in the command is looking for the location of the data. We already set the working directory to C: so the dot is telling R "slow your roll; you don't need to look any further". The second set of quotations in the command is looking for the name of the layer, which in this case is "nhood_2010".

Now we need to prepare our object so that it may be portrayed on a map. R doesn’t know what to do with it in its current form. A few lines of code will transform this caterpillar into a beautiful, map-able butterfly. First, run:

Neighborhoods <- spTransform(Neighborhoods, CRS("+proj=longlat +datum=WGS84"))

spTransform allows us to convert and transform between different mapping projections and datums. This line of code is telling R to convert our Neighborhoods file to longitude/latitude projection and World Geodetic System 1984 datum – a global coordinate (GPS) system used by Google Maps (the initial object was set to a Lambert Conic Conformal projection and a NAD83 datum, as well as a GRS80 ellipsoid). This last bit of information is useful, but you really don’t have to know exactly what it means. Just know that there are a bunch of various coordinate systems that die-hard geography nerds have created (for what I’m sure are good reasons), and all you have to do is smile and remember that we’re essentially just converting our coordinates into a friendly format for integrating with Google Maps (I’m sure I’m going to get heat from one of those geography nerds for diluting it in this way).

Now the fortify command (from the package ggplot2) takes all that wonderful spatial data and converts it into a data frame that R knows understands how to put onto a map.

install.packages('ggplot2')
library(ggplot2)

Neighborhoods <- fortify(Neighborhoods)

Alright meow, we are going to take the map we previously created, BaltimoreMap, and add polygons outlining the neighborhoods from our shapefile. *Side Note: I keep the name of the object the same with each transformation I make. This is a preferential thing. As you are learning, you may wish to name each step differently (e.g. BaltimoreMap1, BaltimoreMap2; Neighborhoods1, Neighborhoods2) so you may go back and look at each one and understand the transformations that take place at each step, which also allows you to identify the area you messed up if you end up receiving an error message along the way.

And now we run the final code:

BaltimoreMap <- BaltimoreMap + geom_polygon(aes(x=long, y=lat, group=group), fill='grey', size=.2,color='green', data=Neighborhoods, alpha=0)
BaltimoreMap

There we have it. She looks good! Notice in the last command we specified that the name of our data is Neighborhoods. This is important. When we set x=long and y=lat, this isn’t just us declaring that we want to use longitude and latitude for our projection; we are telling R that the coordinates for the x (horizontal) and y (vertical) axis of our plot (map) are stored in the columns of our data (Neighborhoods) called ‘long’ and ‘lat’, respectively.

Now play around a bit with the various options for fill, size, color, and alpha (which is the level of transparency from 0 to 1, with the level of opaqueness increasing as you approach 1), as well as the various maptypes and zoom levels from part one. Next session we’ll plot some values (more examples below). Thanks for reading!

All works on this site (spatioanalytics.com) are subject to copyright (all rights reserved) by Paul Bidanset, 2013-2014. All that is published is my own and does not represent my employers or affiliated institutions.