[R] Plotting boundary lines from shapefiles overtop a map of Canada

[Alain Dubreuil]

Hi.  I have a requirement to plot a series of points on a map of Canada along with boundaries defining search and rescue (SAR) regions.  I have been successful in plotting the map of Canada (Lambert projection) and the points, but I have been unable thus far to plot the SAR regions on top of the map.  I'm at the point now where I need help to resolve the issue.

To plot the map of Canada, I have used the following line of code:

      map(database= "worldHires","Canada", ylim=c(39,90), xlim=c(-150,-25), col=alpha("grey90",0.5), fill=TRUE, projection="lambert", param=c(50,65))

Note that the ylim and xlim limits go wider that the actual coordinates of Canada, but that is necessary because the SAR regions go out to sea quite a distance.  Also, I need the map to go all the way to the North Pole.

To plot the points, I have used a "dummy" list of points which I will eventually replace with my real data.  I convert the points to the lambert projection on the fly using the following lines of code:

      lon <- c(-60, -60, -80, -80.1, -90, -105, -140)  #a test longitude vector
      lat <- c(42.5, 42.6, 54.6, 54.4, 75, 68.3, 60)  #a test latitude vector
      coord <- mapproject(lon, lat, proj="lambert", param=c(50,65))  #convert points to projected lat/long 
      points(coord, add=TRUE, pch=20, cex=1.2, col=alpha("red", 0.5))  #plot converted points

As stated, plotting the SAR regions has not worked thus far.  The best I have ever gotten is a square box around the map.  I have data files that list the coordinates of the SAR regions, which is a succession of up to 12100 lat & long points.  A colleague converted those data files into shapefiles defining polygons, with the coordinates already projected to Lambert.  I have tried various options to plot the regions, but none have worked.

Using readOGR:

      region <- readOGR(dsn="C:/myRfolder",layer="mySARshapefile")
      plot(region, add=TRUE, xlim=c(-150,-25),ylim=c(39,90), col=alpha("lightgreen", 0.6), border=TRUE)

Using read.shp and draw.shp:

      region <- read.shp("C:/myRfolder/mySARshapefile.shp")
      draw.shape(shape=region, type="poly", col="red")

Using readShapePoly:

      region <- readShapePoly("C:/ myRfolder/mySARshapefile.shp")
      plot(halRegion, add=TRUE, xlim=c(-150,-25),ylim=c(39,90), col=alpha("lightgreen", 0.6), border=TRUE)

Using readShapeLines after converting the region coordinates to a Lines shapefile instead of a Polygon shapefile:

      region <- readShapeLines("C:/myRfolder/mySARshapefile_lines.shp")
      lines(region, col=alpha("black", 0.6))

I have tried playing with spplot, but I haven't quite understood how this one works yet (gives me an error message: "Error in stack.SpatialPointsDataFrame(as(data, "SpatialPointsDataFrame"),  :   all factors should have identical levels")

I would appreciate any help or insight that you could provide to help me get those boundaries drawn on-top of the country map.

Thanks

Alain Dubreuil
Ottawa, Canada