How to Filter an R Simple Features Collection Using Sf Methods Like St_Intersects()

How to filter an R simple features collection using sf methods like st_intersects()?

Note that st_intersection(, sparse = TRUE) returns a logical matrix, and filter wants a vector. We can get the selection vector by subsetting the matrix:

nc %>%
  filter(st_intersects(., ash_point, sparse = FALSE)[1,])

The . is needed to have nc also be the argument to st_intersects, not only to filter.

It would be nice if the filter.sf method would directly be sensitive to the output of st_intersects, without needing sparse=FALSE and [1,]. I'll put it on some TODO list.

With sf, selecting area features that contain at least one of a selection of point features?

without sample data, this is the best I can come up with.

library( sf )

#find intersecting points/polygons    
intersect <- st_intersection(x = polygons, y = points)

#and go further from there

update

Using the sample data @Spacedman provided in his answer.

library(dplyr)
library(sf)

states %>% 
  #create ID's for the states (if they don't have one already) 
  #state ID should be equal to rownumber (fot the filter later on)
  mutate( id = row_number() ) %>%
  #filter out states that do not have any intersetcions with the points/cities
  filter( id %in% unlist( st_intersects(cities, states) ) ) %>%
  #plot
  mapview::mapview()   

r - Convert output from sf::st_within to vector

Here is how you can get a logical vector from sparse geometry binary predicate:

df$indicator <- st_within(df, box) %>% lengths > 0

or to subset without creating a new variable:

df <- df[st_within(df, box) %>% lengths > 0,]

I cannot test on your large dataset unfortunately but please let me know if it is faster than matrix approach.

How to subset a large SpatialPolygonsDataFrame

If you are using the simple features sf library you can use functions from the tidyverse. Filtering data is a matter of using the filter() function. Notice that you can convert your objects to sf using st_as_sf(). See: https://r-spatial.github.io/sf/reference/st_as_sf.html and How to filter an R simple features collection using sf methods like st_intersects()?

Intersecting Points and Polygons in R

If you do overlay(pts, polys) where pts is a SpatialPointsDataFrame object and polys is a SpatialPolygonsDataFrame object then you get back a vector the same length as the points giving the row of the polygons data frame. So all you then need to do to combine the polygon data onto the points data frame is:

 o = overlay(pts, polys)
 pts@data = cbind(pts@data, polys[o,])

HOWEVER! If any of your points fall outside all your polygons, then overlay returns an NA, which will cause polys[o,] to fail, so either make sure all your points are inside polygons or you'll have to think of another way to assign values for points outside the polygon...

R Overlay points and polygons with a certain degree of tolerance

The example data -

set.seed(1)
library(raster)
library(rgdal)
library(sp)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(0.30*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts1 <- spsample(p2-p, n=3, type="random")
pts2<- spsample(p, n=10, type="random")
pts<-rbind(pts1, pts2)

## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add=TRUE)

Solution using sf and nngeo packages -

library(nngeo)

# Convert to 'sf'
pts = st_as_sf(pts)
p = st_as_sf(p)

# Spatial join
p1 = st_join(pts, p, join = st_nn)
p1

## Simple feature collection with 13 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5.795068 ymin: 49.54622 xmax: 6.518138 ymax: 50.1426
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##   ID_1       NAME_1 ID_2           NAME_2 AREA                  geometry
## 1     1     Diekirch    4          Vianden   76 POINT (6.235953 49.91801)
## 2     1     Diekirch    4          Vianden   76 POINT (6.251893 49.92177)
## 3     1     Diekirch    4          Vianden   76  POINT (6.236712 49.9023)
## 4     1     Diekirch    1         Clervaux  312  POINT (6.090294 50.1426)
## 5     1     Diekirch    5            Wiltz  263  POINT (5.948738 49.8796)
## 6     2 Grevenmacher   12     Grevenmacher  210 POINT (6.302851 49.66278)
## 7     2 Grevenmacher    6       Echternach  188 POINT (6.518138 49.76773)
## 8     3   Luxembourg    9 Esch-sur-Alzette  251 POINT (6.116905 49.56184)
## 9     1     Diekirch    3          Redange  259 POINT (5.932418 49.78505)
## 10    2 Grevenmacher    7           Remich  129 POINT (6.285379 49.54622)

Plot showing which polygons and points are joined -

# Visuzlize join
l = st_connect(pts, p, dist = 1)
plot(st_geometry(p))
plot(st_geometry(pts), add = TRUE)
plot(st_geometry(l), col = "red", lwd = 2, add = TRUE)

EDIT:

# Spatial join with 100 meters threshold
p2 = st_join(pts, p, join = st_nn, maxdist = 100)
p2
## Simple feature collection with 13 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5.795068 ymin: 49.54622 xmax: 6.518138 ymax: 50.1426
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##   ID_1       NAME_1 ID_2           NAME_2 AREA                  geometry
## 1    NA         <NA> <NA>             <NA>   NA POINT (6.235953 49.91801)
## 2    NA         <NA> <NA>             <NA>   NA POINT (6.251893 49.92177)
## 3     1     Diekirch    4          Vianden   76  POINT (6.236712 49.9023)
## 4     1     Diekirch    1         Clervaux  312  POINT (6.090294 50.1426)
## 5     1     Diekirch    5            Wiltz  263  POINT (5.948738 49.8796)
## 6     2 Grevenmacher   12     Grevenmacher  210 POINT (6.302851 49.66278)
## 7     2 Grevenmacher    6       Echternach  188 POINT (6.518138 49.76773)
## 8     3   Luxembourg    9 Esch-sur-Alzette  251 POINT (6.116905 49.56184)
## 9     1     Diekirch    3          Redange  259 POINT (5.932418 49.78505)
## 10    2 Grevenmacher    7           Remich  129 POINT (6.285379 49.54622)

error when using st_intersects, in CPL_geos_binop(st_geometry(x), st_geometry(y), op, par, pattern : Evaluation error: IllegalArgumentException

I checked shapefile nc1 and found features having only 1 point per feature even they're LINESTRING.

Remove these 1-point LINESTRING, the st_intersects operation can be successful.

nc1$cnt2 = stringr::str_count(nc1$geometry, ",")

Here cnt2 is a newly created data.frame column to store number of "," in the geometry column. This can indicate the number of points per feacture.

library('dplyr')
nc3 = filter(nc1, cnt2>1)

We can either do:

res= st_intersects(nc3, nc2)

or:

res = st_join(nc3, nc2, join = st_intersects)

---

Related Topics