3

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?

  • "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." Do I underestand you correctly: yu want to know, where the both datasets overlay and where not? Then I would e.g. use gIntersection from the rgeos package and reclassifiy all the pixels to 1 to get all the overlaying areas. – Iris Nov 10 '15 at 9:50
2

You should show what you have tried in some sort of workflow, ideally with example data. Anyway, you are making things too complicated, I think. I would suggest looking at raster::extract, raster::rasterize to work with the polygons and raster data. If you want to work with the point data derived from the raster data, you can use sp::over or raster::extract

  • Thanks for the help- I think that using 'over' has worked. – Jennifer Taylor Nov 10 '15 at 11:15
2

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. )

1

You should use the function 'over' of the sp package to check if one lat long pair is inside or not a polygon.

1
library(ggplot2)
library(dplyr)
long_and_lat <- species %>% fortify() %>% select(long,lat)
  • 1
    Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format. I think this answer would be more likely to attract upvotes and avoid downvotes if you had provided a few sentences of explanation along with the code. – PolyGeo Feb 1 '17 at 4:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.