R help - Sep 2014 - Plotting boundary lines from shapefiles overtop a map of Canada

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

On 23/09/2014 3:27 a.m., Alain Dubreuil wrote:
> 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)
You don't state which package readOGR() comes from, but it is not part 
of the maps package, which doesn't understand shapefiles, so just using 
plot() on top of a (projected) map() is unlikely to succeed.

I believe what you have to do is go back to your lat/long pairs for your 
SAR regions and use mapproject() to convert them to the coordinates used 
by the plotted projection. Note that you don't need the
proj="lambert"
option when you call mapproject() after a call to map() with a 
projection because the most recent projection (and its parameters) are 
"remembered".  Also I suspect (though untested) is that if you put NA 
pairs in between your lists of projected SAR regions, then you can just 
use lines() to draw them all at once.

Hope this helps,
Ray Brownrigg
> 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
>