[GIS] Extract all the polygon coordinates from a SpatialPolygonsDataframe

coordinatespoint-in-polygonr

I want to do something relatively simply but it seems to be very difficult to do, though I am very new to using R for spatial data so apologies if I've missed something obvious.

I have a list containing 161 shapefiles of bird distributions in the world, they are each in the format of a SpatialPolygonsDataFrame. I also have world climate data in a raster stack format for 8 variables. Both sets of files are in the same projection (lat and long WGS84). I want to create new columns for the climate data with 1's where the species polygon overlaps with the climate data and 0 where it doesn't. I tried various functions to do this, quite a few initial options didn't work, the best option appeared to be:

pnt.in.poly(XY coordinates of climate data points, XY coordinates of species polygon)

I changed my raster of climate data into points to get the XY coordinates. I unified the polygons for each species shapefile and then used:

species@polygons[[1]]@Polygons[[1]]@cords

to get the coordinates for my species polygons as a matrix. However, I discovered that:

species@polygons[[1]]@Polygons[[1]]@coords

wasn't actually getting all the coordinates of the polygons only some of them.

I tried using fortify from the ggplot package, and it gives me a dataframe with coordinates but some of them are for holes in the polygons. I ran that with the pnt.in.poly function excluding the holes but it took an hour to run and when I looked at the output maps after it matched for some of the range but not all of it- there were holes in the distribution where there shouldn't be holes.

I also tried this method that I found somewhere else:

lonlat <- lapply(slot(species, "polygons"), function(x) lapply(slot(x,"Polygons"), function(y) slot(y, "coords")))

which gives me a list with 436 items- but I'm still not sure how to get the coordinates into a matrix form that will work in pnt.in.poly.

I'm sure this should be easier to do! Does anyone know a better way?

Best Answer

There are several ways to get at all of the coordinates for a SpatialPolygonsDataFrame or a SpatialLinesDataFrame.

Each of these methods records the object and the part within the object from the source, albeit with subtle differences and different names.

Say p is a SpatialPolygonsDataFrame. 

#1. coerce to points, then to data frame
# names are Lines.NR Lines.ID Line.NR coords.x1 coords.x2
## start here for polys
lin <- as(p, "SpatialLinesDataFrame")  
## start here for lines
pts <- as.data.frame(as(lin, "SpatialPointsDataFrame"))


#2. convert to matrix 
# names are      object part cump hole         x        y
mpts <- raster::geom(p)

#3. convert to data frame
# names are        long      lat order  hole piece  id group
dpts <- ggplot2::fortify(p)

#4. convert to data frame
## names are   object_ branch_ island_ order_        x_       y_
spts <- spbabel::sptable(p)

Each of these outputs can be used interchangeably as long as you deal with the different name conventions, and for raster::geom remember that it's a matrix. (Actually the coercion one doesn't record the hole/island status so that is less useful depending on what's needed).

You can round trip the last one by

p2 <- spbabel::sp(spts, attr_tab = as.data.frame(p))

(There's an added bonus with spbabel in that you get a 'raster-like' print for the object, but with the data in 'tibble' print form. )

Related Question