Thanks Ray, that works very well! This is my first venture into the R world, and I'm quite impressed with the community's dedication and helpfulness. Thank you again.
Alain -----Original Message----- From: Ray Brownrigg [mailto:ray.brownr...@ecs.vuw.ac.nz] Sent: September-24-14 3:18 PM To: Alain Dubreuil; Ray Brownrigg; William Dunlap Cc: r-help@r-project.org Subject: Re: [R] Plotting boundary lines from shapefiles overtop a map of Canada Hi Alain: The issue is that map() is "too clever" and will only plot those map boundaries that fall within the defined limits, and even the myborder option will not override this. What you can do, with somewhat better results, is something like: map(database= "worldHires", ylim=c(39, 90), xlim=c(-145, -25), col="white", projection="lambert", param=c(50, 65)) map(database= "worldHires", "Canada", ylim=c(39,90), xlim=c(-145,-25), col="grey95", fill=TRUE, projection="lambert", param=c(50, 65), add=T) [Apologies for the HTML formatting in my previous message, I have now turned this off for my home computer.] Hope this helps, Ray On 25/09/2014 12:27 a.m., Alain Dubreuil wrote: > 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. ______________________________________________ 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.
