R: creating a map of selected Canadian provinces and U.S. states

gis, r

Solution

To plot multiple `SpatialPolygons` objects on the same device, one approach is to specify the geographic extent you wish to plot first, and then using `plot(..., add=TRUE)`. This will add to the map only those points that are of interest.

Plotting using a projection, (e.g. a polyconic projection) requires first using the `spTransform()` function in the `rgdal` package to make sure all the layers are in the same projection.

## Specify a geographic extent for the map
## by defining the top-left and bottom-right geographic coordinates
mapExtent <- rbind(c(-156, 80), c(-68, 19))

## Specify the required projection using a proj4 string
## Use http://www.spatialreference.org/ to find the required string
## Polyconic for North America
newProj <- CRS("+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 
            +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

## Project the map extent (first need to specify that it is longlat) 
mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
                  proj4string=CRS("+proj=longlat")),
                  newProj)

## Project other layers
can1Pr <- spTransform(can1, newProj)
us1Pr <- spTransform(us1, newProj) 

## Plot each projected layer, beginning with the projected extent
plot(mapExtentPr, pch=NA)
plot(can1Pr, border="white", col="lightgrey", add=TRUE)
plot(us1Pr, border="white", col="lightgrey", add=TRUE)

Adding other features to the map, such as highlighting jurisdictions of interest, can easily be done using the same approach:

## Highlight provinces and states of interest
theseJurisdictions <- c("British Columbia",
                        "Yukon",
                        "Northwest Territories",
                        "Alberta",
                        "Montana",
                        "Alaska")

plot(can1Pr[can1Pr$NAME_1 %in% theseJurisdictions, ], border="white", 
    col="pink", add=TRUE)

plot(us1Pr[us1Pr$NAME_1 %in% theseJurisdictions, ], border="white", 
    col="pink", add=TRUE)

Here is the result:

Add grid-lines when a projection is used is sufficiently complex that it requires another post, I think. Looks as if @Mark Miller as added it below!

Problem

I am attempting to create a map of selected Canadian provinces/territories and selected U.S. states. So far the nicest maps appear to be those generated with GADM data: http://www.gadm.org/ However, I have not been able to plot the U.S. and Canada on the same map or plot only selected provinces/territories and states. For example, I am interested in Alaska, Yukon, NWT, British Columbia, Alberta, and Montana among others. Also, the U.S. map appears to be split along the international dateline. Can someone please help me to: - plot the aforementioned provinces/territories and states on a single map - avoid having the U.S. split along the International dateline - overlay a latitude-longitude grid - select a specific projection, maybe the polyconic. Maybe spplot does not allow users to specify projections. I did not see an option to select a projection on the spplot help page. I know how to select projections with the map function in the maps package but those maps did not appear to look as nice and I could not plot the desired subset of provinces/territories and states with that function either. I do not know how to begin adding a latitude-longitude grid. However, Section 3.2 of the file 'sp.pdf' seems to address the topic. Below is the code I have come up with so far. I have loaded every map-related package I have stumbled upon and commented out GADM data except for provincial/territorial or state boundaries. Unfortunately, so far I have only managed to plot maps of Canada or the U.S. ``` library(maps) library(mapproj) library(mapdata) library(rgeos) library(maptools) library(sp) library(raster) library(rgdal) # can0<-getData('GADM', country="CAN", level=0) # Canada can1<-getData('GADM', country="CAN", level=1) # provinces # can2<-getData('GADM', country="CAN", level=2) # counties plot(can1) spplot(can1, "NAME_1") # colors the provinces and provides # a color-coded legend for them can1$NAME_1 # returns names of provinces/territories # us0 <- getData('GADM', country="USA", level=0) us1 <- getData('GADM', country="USA", level=1) # us2 <- getData('GADM', country="USA", level=2) plot(us1) # state boundaries split at # the dateline us1$NAME_1 # returns names of the states + DC spplot(us1, "ID_1") spplot(us1, "NAME_1") # color codes states and # provides their names # # Here attempting unsuccessfully to combine U.S. and Canada on one map. # Attempts at selecting given states or provinces have been unsuccessful. # plot(us1,can1) us.can1 <- rbind(us1,can1) ``` Thanks for any help. So far I have made no progress with Steps 2 - 4 above. Perhaps I am asking for too much. Perhaps I should simply switch to ArcGIS and try that software. I have read this StackOverflow post: Can R be used for GIS? EDIT I have now borrowed an electronic copy of 'Applied Spatial Data Analysis with R' Bevand et al. (2008) and downloaded (or located) associated R code and data from the book's website: http://www.asdar-book.org/ I also found some nice-looking GIS-related R code here: https://sites.google.com/site/rodriguezsanchezf/news/usingrasagis If and when I learn how to accomplish the desired objectives I will post solutions here. Although I may eventually move to ArcGIS if I cannot accomplish the objectives in R.

Original source

Related problems