0

I have two very large polygon shapefiles for which I need to calculate the difference between the two shapefiles. One is a shapefile of forest polygons, the other is a shapefile of road polygons. I want to remove the area of the forest polygons that intersects with road polygons, while maintaining each polygon as a separate feature and its data.

After much trial and error, I found that erase using the sp format or using st_erase using the sf format performs the operation that I need.

st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))

However, the processing time is really long. With a small sample of my dataset, I have processing times several minutes long. (And R crashes with a tenth of the dataset). Meanwhile, the function Vector-->Geoprocessing tools--> Difference in QGIS takes 83 seconds with a tenth of the dataset.

erase(forest_sp, roads_sp)
#2.83488 minutes 
forest_sp-roads_sp
#2.8698 minutes
st_erase(forest, roads) 
#58.77sec

Is there a function in R performs this operation quickly? Or is R inferior to QGIS in processing times for this kind of spatial function?

Background

There are multiple questions on gis.stackexchange that deal with this type of problem. Some are with the sp format and suggest gDifference or erase (or its analog polygon1-polygon2):

Erase, gDifference

Clip polygon and retain data?

gDifference

Remove parts of line that fall within a polygon

Identify area within one shapefile not in another shapefile in R

Reverse clipping (erasing ) in R?

Other questions suggest solutions with the sf format: st_difference

Reverse clipping (erasing ) in R?

Using simple features (sf) in R, how do I erase polygons overlapping with another layer

ogr2ogr equivalent of QGIS Union

All state that they perform that same thing (remove intersecting portions) and the official tag definition on gis.stackexchange states that a difference operation removes sections. However, I have found that often they don't. For example, gDifference with byID=FALSE that merges all into one MULTIPOLYGON does remove overlapping sections but gDifference with byID=TRUE that keeps polygons separate includes overlapping sections. With sf, st_difference also keeping overlapping sections and ends up with 10x as many polygons as my original file.

enter image description here

However, erase and the helper function st_erase do work as expected.

enter image description here

Difference

A "difference" operation is a geoprocessing operation that completely removes one overlapping polygon and the intersection of that polygon from the other. It uses one polygon to "take a bite" out of another.

Erase

The action of removing data attributes or features through the processes of positive or negative selection and deletion.

Your Answer

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

Browse other questions tagged or ask your own question.