[R-sig-Geo] Counting points in cells and ploting only cells containing points

Ben Tupper btupper @end|ng |rom b|ge|ow@org
Tue Nov 28 15:37:59 CET 2023


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]]



More information about the R-sig-Geo mailing list