Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: IDs in map2SpatialPolygons

Tags:

r

geospatial

I am trying to use map2SpatialPolygons to create a kernel density plot in a map showing a subset of US States. I keep getting an error saying that "map and IDs differ in length".

I know that this code works (from the help for map2SpatialPolygons):

nor_coast_poly <- map("world", "norway", fill=TRUE, col="transparent", plot=FALSE, ylim=c(58,72))
IDs <- sapply(strsplit(nor_coast_poly$names, ":"), function(x) x[1])
nor_coast_poly_sp <- map2SpatialPolygons(nor_coast_poly, IDs=IDs, proj4string=CRS("+proj=longlat +datum=wgs84"))

this code also works (when mapping the whole US):

usmap <- map('usa', fill=TRUE, col="transparent", resolution=0, plot=FALSE)
uspoly <- map2SpatialPolygons(usmap, IDs=usmap$names, proj4string=CRS("+proj=longlat +datum=WGS84"))

but this code does not:

states.to.plot=c("illinois", "indiana", "ohio")
dmap<-map("state", regions=states.to.plot, col="transparent", plot=FALSE)
dpoly <- map2SpatialPolygons(dmap, IDs=dmap$names, proj4string=CRS("+proj=longlat +datum=WGS84"))

It throws the error:

Error in map2SpatialPolygons(dmap, IDs = dmap$names, proj4string = CRS("+proj=longlat +datum=WGS84")) :
    map and IDs differ in length

How do I get the IDs correct when using map("state"...)?

like image 812
user8026 Avatar asked Jun 19 '12 15:06

user8026


1 Answers

From ?map:

The return value is a list with x, y, range, and names components. ... If fill is FALSE, the x and y vectors are the coordinates of successive polylines, separated by NAs. If fill is TRUE, the x and y vectors have coordinates of successive polygons, again separated by NAs.

So with the default fill = FALSE, you get a bunch of polylines, "a bunch" being more than the three states you wanted to plot, and that's why you get the error above. However, even with the correct number of IDs, you'd still get an error, since map is returning polylines instead of polygons.

With fill = TRUE, you get the polygons that you're after:

dmap<-map("state", regions=states.to.plot, col="transparent", plot=FALSE, fill = TRUE)

dpoly <- map2SpatialPolygons(dmap, IDs = dmap$names,
  proj4string=CRS("+proj=longlat +datum=WGS84"))
like image 89
BenBarnes Avatar answered Sep 18 '22 05:09

BenBarnes