CConley.ca

In addition to the convenience that R offers for data cleaning, analysis and automating reporting, it also has the capacity to complete a variety of mapping (GIS) tasks.  Following are a few R snippets to help get started using the example of plotting schools (as a point file) within their catchment areas (boundaries described by polygons).

Shapefiles: Polygons

The package maptools has a couple of useful functions that will load ESRI shape files into R.  The first is readShapePoly() which, like read.csv, loads polygon .shp files into R as an object.  The second is proj4string() which defines the projection of the shape file:

#load the library

library(maptools)

#load the School_Boundary.shp file into School.bndry object

School.bndry <- readShapePoly("School_Boundary") #loads School_Boundary.shp

#attach the projection used by the School_Boundary.shp file

proj4string(School.bndry) <- "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

This last step, attaching the projection, requires you to:

1. Know the projection that was used in the creation of the shape file
2. Have this projection in proj4 format

If your .shp file also has a.prj file, you can use QGIS to get the proj4text by:

• Opening the shapefile in QGIS
• Right clicking on the boundary file and selecting properties
• Clicking on “Metadata”
• Scrolling to the bottom of the window you will find the proj4 text.  Copy and paste it into R

Working without projection information

If your .shp file does NOT have a .prj file things are a little more challenging.  Here are a few suggestions:

1. If you have another file created by the same organization, check to see what projection it uses.  Organizations tend to be consistent in their file creation and will likely use the same projections from project-to-project and data product-to-data product.
2. Go to the following website, zoom in to your location on the map and click on “PROJ4”: http://www.spatialreference.org/ref/epsg/26912/
3. If you are in Southern Ontario, most of the files I have come across work with UTM NAD 83 zone 17 which is the following in proj4:
   +proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
4. Brute Force: Open a shape file in QGIS that has a projection that is defined.  Open the “mystery” shape file and change it’s projection “on the fly” until you find one that lines up correctly with your first file.  Whichever one lines up is likely the projection you need to use.

Shapefiles: Points

The point file works in a similar manner.  This time the shape file is loaded using the maptools function readShapePoints()

School.point <- readShapePoints("School_points")   #loads School_points.shp

proj5string() is used again to define the projections:

proj4string(School.point) <- "+proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

Dealing with different projections

If the polygon and point file had the same projection I would be ready to prepare for creating a map.  However, the two files have different projections and need to be transformed to have the same projection (whichever one I choose).  In this instance, I’ll create two new objects that transform the files to NAD83 (both are done here for illustrative purposes but the boundary file is already in NAD83):

School.bndry.83 = spTransform(School.bndry,CRS("+init=epsg:26917"))

School.point.83 = spTransform(DDSB.geo.mp,CRS("+init=epsg:26917"))

Plotting the Maps with ggplot2

If you are familiar with the ggplot2 package you will be pleased to know that in addition to plotting histograms, scatterplots etc. it can also plot maps with all the same features. However, to take advantage of ggplot, there is one additional step required with the point file. It needs to be converted to a dataframe:

School.point.83.df <- as.data.frame(School.point.83) #Plotting a Map with ggplot2

With the shape files loaded, projections defined and the point file available as a dataframe the data is now ready to be plotted with ggplot:

 
ggplot(School.bndry,83) +            #Use the school boundary data 
     aes(long,lat, group=group) + 
     geom_polygon() +                #to draw the polygons 
     geom_path(color="white") +      #make the border of the polygons white 
     geom_point(data=DDSB.geo.mp.df, #add the points to the map 
                aes(X, Y, group=NULL, fill=NULL), 
                color="blue", 
                alpha=I(8/10)) #make each point blue with some transparency 

Full Code:

library(ggplot2)

#load the library

library(maptools)

#load the School_Boundary.shp file into School.bndry object
School.bndry <- readShapePoly("School_Boundary")

#attach the projection used by the School_Boundary.shp file
proj4string(School.bndry) <- "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

School.point <- readShapePoints("School_points")   #loads School_points.shp
proj4string(School.point) <- "+proj=utm +zone=17 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

School.bndry.83 = spTransform(School.bndry,CRS("+init=epsg:26917"))
School.point.83 = spTransform(DDSB.geo.mp,CRS("+init=epsg:26917"))

School.point.83.df <- as.data.frame(School.point.83)

ggplot(School.bndry,83) +              #Use the school boundary data
     aes(long,lat, group=group) +
     geom_polygon() +                  #to draw the polygons
     geom_path(color="white") +        #make the border of the polygons white
     geom_point(data=DDSB.geo.mp.df,   #add the points to the map
           aes(X, Y, group=NULL, fill=NULL),
           color="blue", 
           alpha=I(8/10))      #make each point blue with some transparency