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

Alain Dubreuil alain.dubreuil at cae.com
Wed Sep 24 14:27:14 CEST 2014

Bill and Ray,

Thank you to both of you for your input.  I am now able to display the lines.  Both the "polygon" and the "lines" commands that you suggested work.  I do have an extra question for you (or anyone else that can help): when I plot the map of Canada, I specify that I want to show from 25  to 145 degrees W.  However, what shows up is more constrained (probably something like 40 to 125 degrees W).  Is there a way to show just the map of Canada (without the USA) while still showing the desired width?  I assumed that it was constraining the width because of the "Canada" specification in the map command, so I removed it but it ended up showing the entire USA territory (minus Hawaii) in addition to Canada.



-----Original Message-----
From: Ray Brownrigg [mailto:Ray.Brownrigg at ecs.vuw.ac.nz] 
Sent: September-24-14 7:37 AM
To: William Dunlap; Alain Dubreuil
Cc: r-help at r-project.org
Subject: Re: [R] Plotting boundary lines from shapefiles overtop a map of Canada

Thanks Bill for picking this up "while I was sleeping".

An enhancement that I have tested for in the case where the SAR regions are not defined as closed polygons is e.g.:
xValue <- c(105.0, 120.0, 120.0, 105.0, NA, 110, 119, 106) yValue <- c(49.0, 49.0, 60.0, 60.0, NA, 50, 55, 59) polygon(mapproject(y = yValue, x = -xValue))


On 24/09/2014 4:45 a.m., William Dunlap wrote:
> Try:
>    map(database= "worldHires","Canada", ylim=c(39,90), 
> xlim=c(-145,-25), col="grey95", fill=TRUE, projection="lambert",
> param=c(50,65))
>    lines(mapproject(y=yValue, x=-xValue))
> mapproject does care about the sign of the longitude, but if you 
> incompletely reset the projection it messes things up.
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
> On Tue, Sep 23, 2014 at 9:36 AM, William Dunlap <wdunlap at tibco.com> wrote:
>>> testLines <- mapproject(yValue, xValue, proj="lambert", 
>>> param=c(50,65))
>> For starters, if you give the x,y values in reverse order of what the 
>> mapproject function expects you need to label them: y=yValue, 
>> x=xValue.
>> (Also, I would have expected longitudes in the Americas to be 
>> negative, but mapproject doesn't appear to care.)
>> Bill Dunlap
>> TIBCO Software
>> wdunlap tibco.com
>> On Tue, Sep 23, 2014 at 9:17 AM, Alain Dubreuil <alain.dubreuil at cae.com> wrote:
>>> Hi all,
>>> Based on Ray's suggestions, I have tried the following script:
>>> library(mapproj)
>>> library(maps)
>>> resuPdfFileName="C:/linesTest.pdf"
>>> pdf(resuPdfFileName)
>>> # Create a map of Canada in Lambert projection map(database= 
>>> "worldHires","Canada", ylim=c(39,90), xlim=c(-145,-25), 
>>> col=alpha("grey90",0.5), fill=TRUE, projection="lambert", 
>>> param=c(50,65)) # Use test position vectors to draw lines yValue <- 
>>> c(49.0, 49.0, 60.0, 60.0) xValue <- c(105.0, 120.0, 120.0, 105.0) # 
>>> Convert the test vectors into lambert and then draw the lines 
>>> testLines <- mapproject(yValue, xValue, proj="lambert", 
>>> param=c(50,65))
>>> lines(testLines)
>>> dev.off()
>>> The script draws the map of Canada, but fails to draw the lines.  Please let me know what I'm doing wrong because I can't see it.  By the way, not specifying the lambert projection in the call to mapproject yields different results than specifying it which seems contrary to the documentation (?).
>>> Thanks
>>> Alain Dubreuil
>>> Ottawa, Canada
>>> -----Original Message-----
>>> From: Ray Brownrigg [mailto:Ray.Brownrigg at ecs.vuw.ac.nz]
>>> Sent: September-23-14 4:41 AM
>>> To: Alain Dubreuil; r-help at r-project.org
>>> Subject: Re: [R] Plotting boundary lines from shapefiles overtop a 
>>> map of 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
>>> ______________________________________________
>>> 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.

More information about the R-help mailing list