[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.
Thanks
Alain
-----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))
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.
More information about the R-help
mailing list