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.locke using 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 <mga using 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 using r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
[[alternative HTML version deleted]]