Re: [R] Generating maps in R

From: Roger Bivand <Roger.Bivand_at_nhh.no>
Date: Sun, 30 Mar 2008 18:11:35 +0000 (UTC)

Aleksandr Andreev <aleksandr.andreev <at> duke.edu> writes:

>
> Greetings!
>
> I am trying plot some data on a map in R. Here's the scenario.
>
> I have a variable called probworkinghealthy which contains a predicted
> probability of employment for every individual in my sample (about
> 100,000 observations).
> I have another variable, called a001ter, which contains the subject of
> residency in the Russian Federation (akin to a US state) for every
> individual in the sample.
> I have a shape file with the boundaries of all the subjects, called russia.shp.
>
> I can plot boxplots of the probability by Federal Subject using
> plot(probworkinghealthy ~ a001ter). I can also plot the map using
> plot(russia.shp)

So what you need to do is to aggregate the input data frame for individuals, and assign the summary value, here mean, to the correct shapefile geometries. Probably most of what you need is in the maptools package. There is a question about the IDs used to identify the Federal Subject geometries in the shapefile. Assuming that they are named as in the individual data frame, something like:

FS <- readShapePoly("russia.shp", IDvar="a001ter")

will associate the geometries in the SpatialPolygonsDataFrame with the correct ID - change the IDvar argument value to suit the shapefile. If you do not have IDs in the shapefile to match the unique values of a001ter, you will need to create them.

>From there, create a new data frame with row names matching the unique IDs:

agg1 <- tapply(probworkinghealthy, a001ter, mean) agg2 <- data.frame(mean_probworkinghealthy=agg1, row.names=names(agg1))

and check row.names(agg2) and
sapply(slot(FS, "polygons"), function(x) slot(x, "ID"))

for obvious blunders. Merge using:

FS1 <- spCbind(FS, agg2)

(usually takes a number of tries to get the matching correct). The reason for the complications with IDs is that it is easy to merge data with geometries in the wrong order, and having to provide positive matching should help prevent this.

>
> Now, I would like to plot the mean probability of employment (i.e.
> mean(probworkinghealthy)) on a map of Russia using color coding all
> the Federal Subjects. Does anyone know how to do something like that?

Then

spplot(FS1, "mean_probworkinghealthy")

gives a map with default class intervals and colour palette.

I suggest that you also consider assigning a coordinate reference system to the geometries as Russia has a large extent, and the correct interpretation of the geometry coordinates may matter.

Roger

PS. You could consider following up on the R-sig-geo list, or exploring information in the Spatial Task View, or in the list archives.

>
> Much appreciated,
>
> Aleks Andreev
> Duke University
>



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 30 Mar 2008 - 18:39:59 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 30 Mar 2008 - 19:30:25 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive