0

I have two sfc data frames with the same number of rows (individual IDs), but each data frame contains a polygon assigned to the individual. For example, dataframe1 may have a circle assigned neared an individuals' home coordinate, while dataframe2 has a square assigned near their home coordinate.

In some cases, the square and circle overlap for an individual. In some cases, they do not. I want to know, when they do overlap, what is the area of that overlapping space.

At present I've been using st_intersection, but this, of course, returns the area overlapped by all circles (y) for an individual's square (x). I've considered a loop as well, but it seems that st_intersection won't work in cases where x and y do not overlap.

Perhaps a loop with a skip-error function built in? Or a more simpler way...?

  • " I've considered a loop as well, but it seems that st_intersection won't work in cases where x and y do not overlap." how does it not work? It returns an empty geometry. If your code isn't handling that then that's a problem with your code. You are probably very close to a solution with your loop - if you include your code then we might be able to fix it. – Spacedman May 26 at 11:03
0

Make two sf data frames of circles from two sets of 10 points in a unit square:

 set.seed(123)
 p1 = st_cast(st_sfc(st_multipoint(cbind(runif(10),runif(10)))),"POINT")
 p2 = st_cast(st_sfc(st_multipoint(cbind(runif(10),runif(10)))),"POINT")
 b1 =st_buffer(p1, .15)
 b2 =st_buffer(p2, .15)
 b1d = data.frame(id=1:10)
 b1d$geometry = b1
 b1d = st_as_sf(b1d)
 b2d = data.frame(id=1:10)
 b2d$geometry = b2
 b2d = st_as_sf(b2d)

Now compute the intersections by row:

s12 = lapply(1:nrow(b1d), function(i){st_intersection(b1d[i,],b2d[i,])})

ignore the warning messages.

The return list has empty elements if no overlap, geometry elements if there is an overlap:

> s12[[1]]
Simple feature collection with 0 features and 2 fields
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    NA
proj4string:    NA
[1] id       id.1     geometry
<0 rows> (or 0-length row.names)
> s12[[3]]
Simple feature collection with 1 feature and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0.4905068 ymin: 0.5891645 xmax: 0.5589769 ymax: 0.7791115
epsg (SRID):    NA
proj4string:    NA
  id id.1                       geometry
3  3    3 POLYGON ((0.5589769 0.67757...

And if you row-bind them all together you get rid of the empty ones:

> do.call(rbind, s12)
Simple feature collection with 4 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0.401435 ymin: 0.03859731 xmax: 1.013404 ymax: 0.7791115
epsg (SRID):    NA
proj4string:    NA
  id id.1                       geometry
3  3    3 POLYGON ((0.5589769 0.67757...
4  4    4 POLYGON ((0.8638831 0.72135...
5  5    5 POLYGON ((0.8050175 0.03859...
9  9    9 POLYGON ((0.423012 0.250504...

showing you that rows 3,4,5 and 9 in b1d overlap the corresponding rows in b2d. Plotting this object will show you the intersection areas.

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.