[R-sig-Geo] map2SpatialPolygons' 'proj4string' in spmaps

[Roger Bivand]

On Wed, 15 Oct 2008, Sebastian P. Luque wrote:

> On Wed, 15 Oct 2008 19:24:18 +0200 (CEST),
> Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
> [...]
>
>> No, not at all. The coordinates returned by map() if projection= is used
>> are only for plotting, are in arbitrary units, and are only documented in
>> code. map2SpatialPolygons() assumes that the data are in geographical
>> coordinates, and ought perhaps to guess an ellipse and datum based on the
>> age of the dataset (world possibly WGS84, US counties probably NAD27), but
>> allows the user to set the CRS object directly. The projection= argument of
>> map() makes the function call mapproject() with the given projection on the
>> data flowing through, I believe.
>
>> If you do str() of map(projection=) output, you will see the projection
>> name used, but note that the coordinates are arbitrary.
>
> Thanks a lot Roger, so here:
>
> se.baffin <- map("worldHires", xlim=c(-70, -60), ylim=c(62, 69),
>                 resolution=0, interior=FALSE)
>
> or in any of the world* databases, can we safely assume "+proj=longlat
> +ellps=WGS84" (possibly with "+over" in some cases) in the returned
> object?  This would be good info to have in ?map.
>
> One problem I also run into with map2SpatialPolygons() with the above
> map object, is that the ring is not closed.

If the fill=TRUE argument to map() is omitted, the rings are not built:

se.baffin <- map("worldHires", xlim=c(-70, -60), ylim=c(62, 69),
   resolution=0, interior=FALSE, fill=TRUE, plot=FALSE)
baff <- map2SpatialPolygons(se.baffin,
   IDs=se.baffin$names)[match("Canada:Baffin Island", se.baffin$names)]
plot(baff, axes=TRUE)

looks OK to me. Add the CRS later. I don't think that we know definitely 
that WGS84 was used, WDBII seems to have been around since 1985. Try 
overlaying on GE, for example, to get a feel.

> Do you know of any examples showing how one may close such rings, based 
> on a polygon?  Say we have a rectangle over the map object and want to 
> close the ring so we enclose all the land (or water)?

By hand, you need to insert coordinates into the coords slot of the 
Polygon object of interest and re-generate it. A recent question on the 
list about houses east of a road was handled this way (stitch the bounding 
box to the road and use overlay()).

Otherwise look at the unionSpatialPolygons() code for ideas on interfacing 
gpclib classes and methods (also in PBSmapping).

Best wishes,

Roger

>
>
> Cheers,
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no