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.
Thanks Alain -----Original Message----- From: Ray Brownrigg [mailto:ray.brownr...@ecs.vuw.ac.nz] Sent: September-24-14 7:37 AM To: William Dunlap; Alain Dubreuil Cc: r-help@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)) 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 <wdun...@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.dubre...@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.brownr...@ecs.vuw.ac.nz] >>> Sent: September-23-14 4:41 AM >>> To: Alain Dubreuil; r-help@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@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. ______________________________________________ R-help@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.