Help plotting Geographic Data in R using PBSMapping and Shapefiles

Question

Using O'Reilly's Data Mashups in R as inspiration, I'm trying to plot a handful of addresses on a shapefile of Salt Lake County, Utah found here.

I have data frame geoTable:

> geoTable
         address        Y         X EID
1    130 E 300 S 40.76271 -111.8872   1
2    875 E 900 S 40.74992 -111.8660   2
3   2200 S 700 E 40.72298 -111.8714   3
4    702 E 100 S 40.76705 -111.8707   4
5 177 East 200 S 40.76518 -111.8859   5
6    702 3rd ave 40.77264 -111.8683   6
7   2175 S 900 E 40.72372 -111.8652   7
8   803 E 2100 S 40.72556 -111.8680   8

And I've coerced it into an eventData object:

> addressEvents<-as.EventData(geoTable,projection=NA)
> addressEvents
         address        Y         X EID
1    130 E 300 S 40.76271 -111.8872   1
2    875 E 900 S 40.74992 -111.8660   2
3   2200 S 700 E 40.72298 -111.8714   3
4    702 E 100 S 40.76705 -111.8707   4
5 177 East 200 S 40.76518 -111.8859   5
6    702 3rd ave 40.77264 -111.8683   6
7   2175 S 900 E 40.72372 -111.8652   7
8   803 E 2100 S 40.72556 -111.8680   8

So it looks like I've got everything I need to plot-but its not working. When I load the shapefile and plot using

addPoints(addressEvents,col="red",cex=.5)

I'm left looking at an empty shapefile. Additionally, when I try and run findPolys against my eventData object, it returns NULL.

> findPolys(addressEvents,myShapeFile)
NULL

How can I make this work? I was able to complete the O'Reilly tutorial without any problems and am having difficulty figuring out where I'm going wrong here. I dont know if its the shapefile, my data frame, or whateverelse.

Here are the commands I use to import my data and shapefile

slc<-read.table('~/utah.txt',sep=',',header=TRUE,strip.white=TRUE,stringsAsFactors=FALSE)

myShapeFile<-importShapefile("/Users/neil/Downloads/SGID93_DEMOGRAPHIC_CensusTracts2000/SGID93_DEMOGRAPHIC_CensusTracts2000",readDBF=TRUE)

Answer 1

Seems like PBSmapping uses some crude heuristics to work out the projection from the .prj file. (see help(importShapefile)). I personally don't understand all the stuff inside a prj file but using this website www.spatialreference.org I reckon your map matches

http://www.spatialreference.org/ref/epsg/26912/

Whenever I get a new shape file I find it's projection system on this website and then look for the proj4 string, which in this case is "+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

(Like I said I don't know PBSmapping, but you can read this in using maptools as follows)

library(maptools)
sf=readShapeSpatial("SGID93_DEMOGRAPHIC_CensusTracts2000.shp",proj4string=CRS("+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

and then convert to latlong using

library(rgdal)

sftransformed=spTransform(sf,CRS("+proj=longlat"))

and

plot(sftransformed,axes=T)

gives a plot with the right units on the axes.

Not sure if PBSmapping understands a proj4 string? Looks like it doesn't to be honest.

Answer 2

You may also want to look at these related questions, especially at Eduardo's responses:

• plotting-color-map-with-zip-codes-in-r-or-python
• developing-geographic-thematic-maps-with-r