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

Wed Sep 24 13:36:32 CEST 2014

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))

Cheers,
Ray

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.