1

I have a shapefile with points for 6 individuals and another shapefile with 6 polygons (one per individual). Every individual has the same id in both files.

I want to intersect both layers in order to obtain, for each individual, the points within the corresponding polygon. I've tried with only one individual, using the function raster::intersect() and it works perfectly. However, I will have to repeat this process for other animals so I would like to have a loop that could do this automatically (but I am not very familiar with loops).

I've tried several options, and so far this was the only one that worked partially:

for (i in unique(SLO_repr$id)) {
         for (j in unique(HR95_all$id)) {
                 if ((!is.na(SLO_repr$id[i])) == (!is.na(HR95_all$id[j]))) {
                         a <- raster::intersect(SLO_repr, HR95_all)
                         writeOGR(a, dsn= getwd(), paste0("points_withinHR", "_", i), driver = "ESRI Shapefile",overwrite_layer=TRUE)}
         }
 }

The problem with this one is that, in each shapefile, all points within each polygon are included. So, the raster::intersect() function worked, but it seems I didn't select the i and j properly...?

To build this loop, I followed this question: Loop to check multiple polygons overlap in r

I was also wondering if I really need to have a nested loop? Because all id's are the same, so in the condition, when it gets on i for SLO_repr$id, it should select the same value for HR95_all$id?

So I've tried this loop:

 for (i in unique(SLO_repr$id)) {
                 if ((!is.na(SLO_repr$id[i])) == (!is.na(HR95_all$id[i]))) {
                         a <- raster::intersect(SLO_repr[[i]], HR95_all)
                         writeOGR(a, dsn= getwd(), paste0("points_withinHR", "_", i), driver = "ESRI Shapefile",overwrite_layer=TRUE)}
         }

But I get this error: "Error in as.vector(y) : no method for coercing this S4 class to a vector"

Any ideas about this?

1 Answer 1

2

You can just use one loop, where you subset the point layer by unique id. Then intersect these subsetted points by polygon with same id. Code below suggesting that there is one polygon for each unique point id, otherwise subset it in same way as points.

for (id in unique(points$id)){
    id_points <- points[points$id==id,]
    sp_points <- intersect(id_points,polygons[id,])
    writeOGR(sp_points,dsn=getwd(), layer=paste0("points_withinHR_", id), driver = "ESRI Shapefile",overwrite_layer=TRUE)
}
1
  • This worked perfectly, thank you!!!
    – mto23
    May 1, 2020 at 23:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.