Oh, that's an interesting approach; I haven't seen `count()` used like that - nice!.
Here's an alternative that uses stars as a raster container, and the `st_cells()` function to determine which cell each point belongs to. ### START suppressPackageStartupMessages({ library(stars) library(sf) library(dplyr) }) #' Given cell indices, convert to col-row addresses #' @export #' @param x stars object #' @param cells num, cell addresses #' @return matrix of [col, row] colrow_from_cells <- function(x, cells){ d = dim(x) if (any(cells > prod(d))) stop("cell indices must <= nrow(x) * ncol(x)") col = ((cells-1) %% d[1]) + 1 row = floor((cells - 1) / d[1]) + 1 cbind(col, row) } R = system.file("tif/olinda_dem_utm25s.tif", package = "stars") |> stars::read_stars() set.seed(1345) pts = sf::st_bbox(R) |> sf::st_as_sfc() |> sf::st_sample(10000) |> sf::st_as_sf() |> sf::st_set_geometry("geometry") index = stars::st_cells(R, pts) xy = colrow_from_cells(R, index) # augment the pts with index, row and column # group by index, get the count per group pts = dplyr::mutate(pts, index = index, row = xy[,1], col = xy[,2]) |> sf::st_drop_geometry() |> dplyr::group_by(index) |> dplyr::group_map( function(tbl, key){ dplyr::slice(tbl, 1) |> dplyr::mutate(n = nrow(tbl)) }, .keep = TRUE) |> dplyr::bind_rows() # make a "blank" copy of the raster, assign the counts to each cell C = R*0 names(C) = "point count" C[[1]][pts$index] <- pts$n plot(C, axes = TRUE) ### END On Wed, Nov 22, 2023 at 7:00 AM Dexter Locke <dexter.lo...@gmail.com> wrote: > Here is one approach. Kind of a silly example with just one grid cell, but > would work with more polygons and more points. > > > library(sf) > library(tidyverse) > > # make some points data > # modified example from > # https://r-spatial.github.io/sf/reference/geos_binary_pred.html > (pts <- > st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, > 2.5))) |> > st_as_sf()|> > rowid_to_column(var = 'pts_id')) > > (pol <- > st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0)))) |> > st_sfc() |> > st_as_sf()|> > rowid_to_column(var = 'pol_id') > ) > > # look at the data, crude 'map' > plot(pts); plot(pol, add = TRUE) # take a look > > # perform the intersection > pts_pol_int <- > pts |> > st_intersection(pol) # only going to retain the overlaps (drops > non-intersecting elements) > > # count the overlaps > cont_pts_pol_int <- > pts_pol_int |> > st_drop_geometry() |> # pulls out the data frame (or tibble) > group_by(pol_id) |> > count() > > cont_pts_pol_int # could be joined back to pol > > HTH, DHL > > > On Wed, Nov 22, 2023 at 6:39 AM MIGUEL GOMEZ DE ANTONIO <m...@ccee.ucm.es> > wrote: > > > Dear community, > > > > I want to plot only the cells of a grid that contains points. > > > > I have a set of points (of class SpatialPointsDataFrame) and a grid (of > > class SpatialPolygonDataFrame). > > > > I am interested in ploting the set of cells of the grid that contains > > points: > > > > 1. Count how many points fall in each cell of the grid. > > > > 2. Plot only the cells that contains points. > > > > I tried: > > > > library(sf) > > > > Points.sf=st_as_sf(points) > > > > Grid.sf=st_as_sf(grid) > > > > A=intersects(grid.sf, points.sf) > > > > apply(A,1,any) > > > > There I get the cells that contains points but: > > > > 1. How can I count how many points fall in each cell > > > > 2. How could I plot the set of cells that contains points? > > > > > > Thank you very much for your help, > > > > > > Miguel Gómez de Antonio > > Profesor Titular de Universidad > > Dpto. Economía Aplicada, Pública y Política > > Universidad Complutense de Madrid > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo