16 Sept 2013

R GIS: Polygon Intersection with gIntersection{rgeos}

A short tutorial on doing intersections in R GIS. gIntersection{rgeos} will pick the polygons of the first submitted polygon contained within the second poylgon - this is done without cutting the polygon's edges which cross the clip source polygon. For the function that I use to download the example data, url_shp_to_spdf() please see HERE.


library(rgeos)
library(dismo)

URLs <- c("http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_gew_pl.zip",               # all water bodies in Tyrol
          "http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_tseepeicher_pl.zip")       # only artificial..

y <- lapply(URLs, url_shp_to_spdf)
z <- unlist(unlist(y))
a <- getData('GADM', country = "AT", level = 2)

b <- a[a$NAME_2=="Innsbruck Land", ]                                               # political district's boundaries
c <- spTransform(b, z[[1]]@proj4string)                                            # (a ring polygon)    
z1_c <- gIntersection(z[[1]], c, byid = TRUE)                                      
z2_c <- gIntersection(z[[2]], c, byid = TRUE)

plot(c)
plot(z1_c, lwd = 5, border = "red", add = T)
plot(z2_c, lwd = 5, border = "green", add = T)
plot(z[[1]], border = "blue", add = T)              # I plot this on top, so it will be easier to identify
plot(z[[2]], border = "brown", add = T)

4 comments :

  1. Nice blog, a fantastic recycle-bin!!! I'll link your blog in my notasr blog!

    ReplyDelete
  2. Works fine, but with this method, you loose attached informations in attribute table ...

    Is there a method to have a similar result of an intersect in ArcGIS or QGIS (intersect geometries + corresponding attributes of initial data) ??

    ReplyDelete
    Replies
    1. This answer looks like it would join attributes into the SpatialPolygons created by gIntersection, but my attempt didn't quite match it back and haven't debugged it yet:
      http://gis.stackexchange.com/questions/90575/join-grid-or-polygon-layers-in-r

      Delete