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

[William Dunlap]

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