[R-sig-Geo] map to shapefile with projection info

Roger Bivand Roger.Bivand at nhh.no
Wed Dec 6 20:28:42 CET 2006

On Wed, 6 Dec 2006, Sebastian P. Luque wrote:

> Hi,
> Trying to extract a map (package maps) to a shapefile, IIUC the procedure
> is to get it into a SpatialPolygonsDataFrame and then use
> writePolyShape().  This seems to require an intermediate step through
> SpatialPolygons.  However, I'm running into some problems:
> R> canada <- map("world2Hires", "Canada", plot=FALSE)
> R> canada_sp <- map2SpatialPolygons(canada, IDs=canada$names,
> +                                  proj4string=CRS("+proj=longlat +ellps=wgs84"))
> Error in map2SpatialPolygons(canada, IDs = canada$names, proj4string = CRS("+proj=longlat +ellps=wgs84")) : 
> 	map and IDs differ in length

canada <- map("world2Hires", "Canada", fill=TRUE, plot=FALSE)

fill=TRUE retrives polygons, otherwise you get line segments. For this 
example it now works. You can also make several retrieved island belong to 
the same Polygons object by giving them the same ID. There is an example 
on the pal2SpatialPolygons help page.

> How can these two components be lined up?
> I would also appreciate some suggestions on how to write a projection file
> for the shapefile.  I can see this has been discussed in the past, but at
> that time there was no solution.

In rgdal (the PROJ.4 library is needed):

names <- sapply(slot(canada_sp, "polygons"), function(i) slot(i, "ID"))
canada_spdf <- SpatialPolygonsDataFrame(canada_sp, 
  data=data.frame(names=names, row.names=names))
writeOGR(canada_spdf, dsn=".", layer="canada", driver="ESRI Shapefile")

where canada.prj is:



> 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

More information about the R-sig-Geo mailing list