Map in R -- Examples, part 2

Maps in R – Examples (Part 2)

This part of the examples web page describes map projection done using the rgdal package (and the sp class and maptools packages).

3. Projections

The sp class and maptools package provide a mechanism for doing projected maps. (Note that the projection parameters used in the example here are not really appropriate for the area being plotted, but were chosen to make the fact that the data are projected evident.)

First, load the rgdal package

library(rgdal)

This example plots a projected map of Oregon climate station data where the data are contained in the in the orstations shapefile. The first block of code does some set up (class intervals and colors, etc.)

# equal-frequency class intervals -- spplot & projected
plotvar <- orstations.shp@data$tann # gets data from shapefile .dbf
nclr <- 8
plotclr <- brewer.pal(nclr,"PuOr")
plotclr <- plotclr[nclr:1] # reorder colors
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
basemap <- list("sp.lines", orotl.shp, fill=NA)

The next block adds “projection string” (i.e. the name of the projection that the data are currently in to the shape file (in case it’s not already there), and does the projection. First a long string of information in specified in the first statement, then the spTransform() function does the actual conversion of latitudes and longitudes into x, y coordinates. Next, the gridlines() function creates a graticule in latitude and longitude, and the spTransform() function projects those. Finally, the spplot() function makes the map.

Here’s the map:

Note that the map is a little different from that made using the standard plot() function. spplot()uses the “lattice package” for plotting, and (in these examples) automatically adds a legend.

Here’s a second example, Oregon county census data, where the data being plotted comes from the orcounty.shp shape file

# project Oregon county data
aea.proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-110
    +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m"
orcounty.shp.proj <- spTransform(orcounty.shp, CRS(aea.proj))

# equal-frequency class intervals -- spplot & projected
plotvar <- orcounty.shp@data$AREA
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu")
class <- classIntervals(plotvar, nclr, style="equal")
colcode <- findColours(class, plotclr)
spplot(orcounty.shp.proj, "AREA",
    col.regions=plotclr, at=round(class$brks, digits=1))

# equal-width class intervals -- spplot & projected
plotvar <- orcounty.shp@data$AREA
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu")
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)

spplot(orcounty.shp.proj, "AREA",
    col.regions=plotclr, at=round(class$brks, digits=1))

Note that the code for producing maps using the spplot() function is much like that using plot(): some setup work first, then determination of the actual colors or symbols that will be plotted, then finally creation of the map.

4. Examples using the maps package

The R maps package provides a means of mapping data that are not necessarily components of a shapefile. The package provides a way of plotting choropleth maps using polygons that it contains (U.S. states and counties, countries of the world), and can use it’s internal polygons to provide unfilled basemaps for point data. The following example replicate two maps done using the maptools package.

Oregon county census data -- bubble plots of orcounty.csv data

library(maps)

# look at the orcounty.csv data
attach(orcountyp)
summary(orcountyp)

# check names for proper order
map("county","oregon", names=T, plot=F)
print(orcountyp$Name)

# bubble plot equal-frequency class intervals
plotvar <- orcountyp$Area
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu")
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
max.symbol.size=10
min.symbol.size=1
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
symbol.size <- ((plotvar-min(plotvar))/
(max(plotvar)-min(plotvar))*(max.symbol.size-min.symbol.size)
+min.symbol.size)

map("county", "oregon", xlim=c(-125,-114), ylim=c(42,47))
map.axes()
points(orcountyp$Longitude, orcountyp$Latitude, pch=16, col=colcode,
cex=symbol.size)
points(orcountyp$Longitude, orcountyp$Latitude, cex=symbol.size)
text(-120, 46.5, "Equal-Frequency Class Intervals")
legend(locator(1), legend=names(attr(colcode, "table")),
fill=attr(colcode, "palette"), cex=0.6, bty="n")

The map() function does the main work of making the map, but has some other applications. When the data are contained in a shape file, the correct attributes (one hopes) are automatically associated with each point or polygon. In this example, the data come from a .csv file, and so the first thing to do is to make sure that the right lines in the .csv file are being matched with the right polygons that the map() function uses. The two lines of code

map("county","oregon", names=T, plot=F)
print(orcountyp$Name)

print out the names of the polygons, and the county names in the spreadsheet. The map() function assumes that the data/polygons are in alphabetical order. Next the usual stepup work is done, the basemap is plotted using the map() function, and axes are labeled using the map.axes() function. Finally, the “bubbles” are added using the points() function. Note that this map is unprojected.

5. Basemaps using the maps package

The maps package provides a means of constructing basemaps for plotting the locations of points, which can be decorated with text, symbols, and so on -- most of the things that be done on scatter plots.

The following example plots the location of large cities. A built-in data set world.cities is made available using the data() function, and some information describing each of the cities in the cities dataframe is gotten using the match() function, with looks at the name of a city in the cities dataframe, and looks in world.cities for the information (latitude and longitude in this case).

# map of large cities

data(world.cities) # make the world cities location data set from the maps package available

# match the large cities with those in the database

m <- match(paste(tolower(as.character(cities$City)),tolower(as.character(cities$Country))),
paste(tolower(world.cities$name),tolower(world.cities$country.etc)))

# assign the world.cities location information to the large cities

big.cities <- NULL

big.cities$name <- cities$City

big.cities$long <- world.cities$long[m]

big.cities$lat <- world.cities$lat[m]

big.cities

# plot the map

map("world")

map.axes()

points(big.cities$long,big.cities$lat, col="blue")

text(big.cities$long, big.cities$lat, big.cities$name, col="red", cex=.5)

The next example generates a projected map of large cities. The projection is done using the mapproject() function.

# map of large cities
# match the large cities with those in the database

m <- match(paste(tolower(as.character(cities$City)),tolower(as.character(cities$Country))),
paste(tolower(world.cities$name),tolower(world.cities$country.etc)))

big.cities <- NULL
big.cities$name <- cities$City
big.cities$long <- world.cities$long[m]
big.cities$lat <- world.cities$lat[m]

# map projection information
proj.type <- "azequalarea"
proj.orient <- c(90,0,30)

# plot the map
map("world", proj=proj.type, orient=proj.orient, resolution=0, wrap=T)
map.grid(col="black", labels=F, lty=1)
proj.coords <- mapproject(big.cities$long,big.cities$lat, proj=proj.type, orient=proj.orient)
points(proj.coords, col="blue")
text(proj.coords, labels=big.cities$name, col="red", cex=0.5)

6. Further examples illustrating map projection using the maps package

The maps package can provide projected base maps, which can provide less distorted views of a data set. Here are some examples that plot the locations of the Oregon climate-station data

unprojected maps (left map below)

# unprojected
map("county", "oregon", fill=F)
points(orstationc$lon, orstationc$lat)
text(orstationc$lon, orstationc$lat, labels=orstationc$station, col="red",
cex=.8)

# unprojected, add surrounding states
map("county", c("oregon","washington"), xlim=c(-125,-116), fill=F)
map("state", "oregon", fill=F, col="grey", lwd=3, add=T)
points(orstationc$lon, orstationc$lat)
text(orstationc$lon, orstationc$lat, labels=orstationc$station, col="red",
cex=.8)

projected maps (right map below)

# projected
proj.type <- "albers"
proj.stdlats <- c(29.5, 45.5)
proj.orient <- c(90,-120,0)
map("county", c("oregon","washington"), proj=proj.type, par=proj.stdlats, orient=proj.orient, fill=F)
orstationc.xy <- mapproject(orstationc$lon, orstationc$lat, proj=proj.type,
orient=proj.orient, par=proj.stdlats)
map("state", "oregon", proj=proj.type, par=proj.stdlats,
orient=proj.orient, fill=F, col="grey", lwd=3, add=T)
points(orstationc.xy)
text(orstationc.xy, labels=orstationc$station, col="red", cex=.8)

7. Basemap shapefile generation using the maps package

The map2SpatialLines() function in the maptools function can be used to transform lines extracted from the maps package into sp() package-compatible format to provide basemaps for plotting other data sets. Here is an example for the Pacific Northwest.

#library(maptools) # also loads sp package
#library(maps)

# extract county outlines from maps() database
pnw.outlines <- map("county", c("oregon","washington", "california",
    "nevada", "idaho"),
xlim=c(-124.5, -116.0), ylim=c(41.0, 47.0))

# prune the lines to Washington, Oregon, and Northern California extent
pnw.outlines <- pruneMap(pnw.outlines, xlim=c(-125.0, -115.0), ylim=c(41.0,
     47.0))

# convert to sp lines
pnw.outlines.sp <- map2SpatialLines(pnw.outlines, proj4string=
     CRS("+proj=longlat"))

plot(pnw.outlines.sp, col="gray", lwd=2)
degAxis(1, at=seq(-125.,-116., by=1.))
degAxis(2, at=seq(42.,47., by=1.))

# project the outlines
aea.proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5
     +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m"
pnw.outlines.proj <- spTransform(pnw.outlines.sp, CRS(aea.proj))

# generate and project gridlines
grat <- gridlines(pnw.outlines.sp, easts=seq(-127,-115,by=1),
     norths=seq(40,47,by=1), ndiscr = 20)
grat.proj <- spTransform(grat, CRS(aea.proj))

plot(pnw.outlines.proj, col="gray", lwd=2)
lines(grat.proj, col="blue", lty=3)

Here’s what you get:

There’s a lot more that can be accomplished, but remember that R is not really a cartographic design-type of package (if there even is one). The main utility of the mapping tools is there ability to visualize data rapidly.