1

I downloaded two shapefiles from IUCN (available here and here) to determine the range overlap between two species. When I bring the shapefiles into R, I am able to plot them on a map together, however, when I try to calculate percent overlap, I get this error:

 p <- readOGR(getwd(), layer="species_7928")
 w <- readOGR(getwd(), layer="species_21314")
 EPFU<- gArea(SpatialPolygons(p@polygons))
 intersected.area.tabr <- gArea(gIntersection(SpatialPolygons(w@polygons), 
 SpatialPolygons(p@polygons)))
 (intersected.area.tabr/EPFU)*100

  Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id,drop_lower_td,
  unaryUnion_if_byid_false,  : TopologyException: Input geom 0 is
  invalid: Ring Self-intersection at or near point -82.517899999999997 
  8.2930355900000006 at -82.517899999999997 8.2930355900000006

I've tried the suggested workarounds with gBuffer (width=0) and gUnaryUnion, but these also return errors.
I was also able to take a separate map provided by a friend from a different project and each of these files individually worked when I tried to calculate the intersection. I can open these files together in ArcMap. I've run the check geometry tool there and the output table was empty. Is this occurring because there are multiple elements (discontinuous portions of the map)? I do not necessarily have to do this in R, if there is a straightforward way to do the same in ArcMap, but this is my first experience with any sort of GIS data.

4

Try the new simple features package (sf):

library("sf")
library("magrittr")

p <- st_read("species_7928.shp") %>%
     st_transform("+init=epsg:2163")  # project geometry
                                      # e.g. US National Atlas Equal Area
w <- st_read("species_21314.shp") %>%
     st_transform("+init=epsg:2163")
res <- st_intersection(st_geometry(p),
                       st_geometry(w))

plot(st_geometry(p), axes=TRUE)
plot(st_geometry(w), col="lightblue", add=TRUE)
plot(res, col="#FF000055", add=TRUE)

plot

sum(st_area(p))
# 1.279e+13 m^2
sum(st_area(w))
# 1.386e+13 m^2
sum(st_area(res))
# 6.108e+12 m^2
  • This worked perfectly. I can divide the total area of the intersection by the species area to get percent overlap. Still, I would like to know what it is about the files that make them give these errors when others do not. – user89092 Jan 4 '17 at 21:24
2

You have invalid (self intersecting) geometries. You can see that by running gIsValid on your SpatialPolygons* objects. They need to be fixed (the zero buffer approach does the trick). Also, for area calculations to be meaningful, SpatialPolygons* should be projected, preferably using an equal area projection.

See the approach below using sp (imported via raster)

library(raster) # For convenience, shapefile function and show method for Spatial objects
library(rgeos)
library(magrittr)

# Lambert equal area with origin on the center of your area
prj <- '+proj=laea +lat_0=10 +lon_0=-81 +ellps=WGS84 +units=m +no_defs'

sp1 <- shapefile('species_21314/species_21314.shp') %>%
      spTransform(CRS(prj)) %>%
      gBuffer(byid=TRUE, width=0)

sp2 <- shapefile('species_7928/species_7928.shp') %>%
      spTransform(CRS(prj)) %>%
      gBuffer(byid=TRUE, width=0)

sp3 <- gIntersection(sp1, sp2)

gArea(sp3)

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.