4

I have two polygon vector layers one of which has more than 100,000 polygons. There are areas of overlap in the two vector layers however they have different attributes. In the final layer I want the attributes to come across from both the polygon layers especially where the polygons in the two layers overlap. I have to do further processing on this attribute intersect information hence it is critical.

I already know how to do this in R (using the union function as part of the raster package) because that's where I'm carrying out most of the post-processing of this data however it is extremely slow. Slow to the point that it continues to do the union even after 2 hours and which it crashes. In QGIS, it takes less than a couple of minutes but I want to run it as a routine process so want to avoid using QGIS.

I have since then tried to do this using ogr2ogr and I am not getting the desired result. For example, I have followed this particular thread and the result is not what I am after. What I am getting at the end is a massive shapefile with the attributes of the two layers combined into 1 but no mapping of the attributes where the intersection has occurred as such. It almost looks like the two shapefiles have been merged without any intersection occurring between them which is what I am after.

Isn't there a union function in ogr2ogr similar to the one in QGIS? Or can anyone suggest how I can expedite this union process in R?

4

If you are interested in a "Spatial OR" between the layers 1 and 2 (i.e. intersection + difference 1-2 + difference 2-1), the union of the three elementary operations can be made with ogr2ogr in this way:

ogr2ogr union.shp layer1.shp -dialect SQLite -sql "select ST_Intersection(a.geometry, b.geometry) as geometry, a.*, b.* from layer1 a, 'layer2.shp'.layer2 b where ST_Intersects(a.geometry, b.geometry) union select ST_Difference(a.geometry, b.geometry) as geometry, a.*, b.* from layer1 a, 'layer2.shp'.layer2 b union select ST_Difference(b.geometry, a.geometry) as geometry, a.*, b.* from layer1 a, 'layer2.shp'.layer2 b"

Although, I suggest to use directly a SpatiaLite or PostGIS db with a spatial index in order to improve the speed of the whole task.

  • 1
    Can't you use some ST_Union command ? – gisnside Aug 10 '17 at 7:49
  • 1
    @gisnside Unfortunately ST_Union dissolves all the feature types (intersections and differences) with the same attributes. Instead, splitting in three elementary geoprocesses allows to isolate each different feature type. This solution should be further refined in order to have null values where there aren't overlaps. – Antonio Falciano Aug 10 '17 at 8:51
  • @afalciano- Thanks a lot for your answer and so sorry for my late response. I tried your solution and it was still slow. So far, the fastest way to get what I am after is using the RQGIS plugin from within R. – VGu Sep 13 '17 at 23:48
0

I tried a solution where I separated the union of the two polygon shapefiles as the combination of an intersection and a difference. On my little example repeated 100 times, this seems 3x faster than the raster::union. Below is the code for my little example, would you try it and tell me if this is faster for your datasets.

library(sp)
library(rgdal)
library(raster)

# Create Polygons
p1 <- rbind(c(-180,-20), c(-140,55), c(-50, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
data1 <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(p1)), 1), Polygons(list(Polygon(p2)), 2))),
                                  data = data.frame(col1 = 1:2), match.ID = FALSE)

p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
data2 <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(p3)), 3))), 
                                  data = data.frame(col2 = "test"), match.ID = FALSE)

Here we retrieved the combined geometry. Disaggregation allows to get all subpolygons separately.

# Geometry
inter.sp <- rgeos::gIntersection(data1,data2) %>% sp::disaggregate()
diff.sp <- rgeos::gSymdifference(data1,data2) %>% sp::disaggregate()
# Union of geometries
comb.sp <- bind(inter.sp, diff.sp)

This only allows to retrieve the geometry and not data that interest you, so I used function sp::over to get it. The previous disaggregation has full sense here.

# Get data
comb.sp.df <- SpatialPolygonsDataFrame(comb.sp,
  data = cbind(over(comb.sp, data1), over(comb.sp, data2)))

The output is the targeted SpatialPolygonsDataFrame

>comb.sp.df
class       : SpatialPolygonsDataFrame 
features    : 5 
extent      : -180, 160, -60, 60  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 2
names       : col1, col2 
min values  :    1, test 
max values  :    2, test 

Note: The same trick with library sf is not quicker, moreover, there seems to be a problem with st_difference that I will notify as a bug.

0

Here's a potential solution using R package sf. If I understand you correctly you want the complete set of intersection between the two polygon sets plus their non-overlapping parts.

library(sf)
library(mapview) # for the data and the overlap checking

# pol1 = st_read("/datadisk/bu_lenovo/software/testing/mapview/switzerland/landuse.shp")
pol1 = franconia
# pol2 = pol1[1:200, ]
pol2 = pol1
st_geometry(pol2) = st_geometry(pol2) + 0.5 # offset pol2 geometry by half a degree
st_crs(pol2) = st_crs(pol1) # the above loses crs information so reset

# do they overlap?
viewExtent(pol1) + viewExtent(pol2)

st_or = function(x, y, dim = 2) {

  # st_erase to get the remainder of the intersection (taken from ?st_difference)
  st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))

  # we need st_dump to extract polygons from a potential GEOMETRYCOLLECTION
  st_dump = function(x, dim) {
    dims = sapply(x, st_dimension)
    x[dims == dim, ]
  }

  # get overlap via intersection
  overlap = sf::st_intersection(x, y)

  # extract polygons (if dimm = 2)
  overlap = overlap[st_dimension(overlap) == dim, ]
  gc = which(sapply(seq(nrow(overlap)), function(i) {
    inherits(overlap[i, ], "GEOMETRYCOLLECTION")
  }))
  if (length(gc) > 0) {
    dmp = st_dump(overlap, dim = dim)
    overlap = rbind(overlap[-gc, ], dmp)
  }

  # get the non-intersecting parts and set missing attribute to NA
  diff = rbind(st_erase(x, y), st_erase(y, x))
  diff[, setdiff(names(overlap), names(diff))] = NA

  # return combined geometry set with attributes
  return(rbind(overlap, diff))

}

system.time({
  combo = st_or(pol1, pol2)
})

viewExtent(combo) + viewExtent(pol1) + viewExtent(pol2)

The commented out pol1 variant (and the corresponding commented out pol2 call) uses the OSM landuse layer of Switzerland (about 97000 Polygons) downloaded from here - the .shp.zip file. Processing the complete set (pol1) with a subset of 200 Polygons of the same set (pol2) takes about 5 min on my machine (and harly uses any RAM), though it is really hard to infer any sort of performance baseline from this as it highly depends on the complexity of the features.

I have done a quick profiling of the st_or function and I think the bottleneck is st_erase. I think there is room for improvement, but I'd suggest to try this implementation first to see whether it produces the result you are after.

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.