R – Count Multilinestring Elements Crossing a Raster Cell

countrrasterrasterization

I have a shapefile with several MULTILINESTRING geometries, representing the borders of some species' distributions, and I want to count the number the number of species occurring in each cell.

I used the rasterize function from the terra package:

count.sp <- terra::rasterize(sp.lines, rast.w, fun = "count")

As an example, the lines in the figure below represent the borders of the distribution of a given species "X":

In the example, species' distribution was subdivided into several disjoint parts crossing the same raster pixel. With the code above, I got as a result the number of parts in each pixel (8 and 5, respectively), whereas I would like to get all the parts belonging to the same species and crossing the same pixel be counted as a single occurrence (i.e., I would like to get 1 in both pixels in the example below).

I also tried the function "sum" (with no field, and I obtained 208 and 130, I don't understand what these values represent)

You can find some example data here (1 shapefile with just 2 species and the raster)

What I am missing?

Best Answer

Following @JeffreyEvans suggestion, I achieved it using an "extract" approach:

First, I assigned an unique ID to each raster pixel:

values(rast.w) <- 1:ncell(rast.w)

Then, I used the extract function from the package terra to get the ID of the cells crossed by each species:

sp <- terra::extract(rast.w, sp.lines)

Then, I converted the resulting list into a dataframe with the ID of each cell (first column) and the number of species occurring in each of them (second column, aggregated counting)

count <- as.data.frame(table(unlist(sp)))

Finally, I created a new raster from the first one, replacing ID values with the counting result:

newraster <- rast.w
newraster <- subs(rast.w, count, by = "test")

Related Solutions

Raster Clumping – Creating Randomly Shaped Clumps of Cells from Seed Pixels

I will offer an R solution that is coded in a slightly non-R way to illustrate how it might be approached on other platforms.

The concern in R (as well as some other platforms, especially those that favor a functional programming style) is that constantly updating a large array can be very expensive. Instead, then, this algorithm maintains its own private data structure in which (a) all the cells that have been filled so far are listed and (b) all the cells that are available to be chosen (around the perimeter of the filled cells) are listed. Although manipulating this data structure is less efficient than directly indexing into an array, by keeping the modified data to a small size, it will likely take much less computation time. (No effort has been made to optimize it for R, either. Pre-allocation of the state vectors ought to save some execution time, if you prefer to keep working within R.)

The code is commented and should be straightforward to read. To make the algorithm as complete as possible, it makes no use of any add-ons except at the end to plot the result. The only tricky part is that for efficiency and simplicity it prefers to index into the 2D grids by using 1D indexes. A conversion happens in the neighbors function, which needs the 2D indexing in order to figure out what the accessible neighbors of a cell might be and then converts them into the 1D index. This conversion is standard, so I won't comment on it further except to point out that in other GIS platforms you might want to reverse the roles of column and row indexes. (In R, row indexes change before the column indexes do.)

To illustrate, this code takes a grid x representing land and a river-like feature of inaccessible points, starts at a specific location (5, 21) in that grid (near the lower bend of the river) and expands it randomly to cover 250 points. Total timing is 0.03 seconds. (When the size of the array is increased by a factor of 10,000 to 3000 rows by 5000 columns, the timing goes up only to 0.09 seconds--a factor of only 3 or so--demonstrating the scalability of this algorithm.) Instead of just outputting a grid of 0's, 1's, and 2's, it outputs the sequence with which the new cells were allocated. On the figure the earliest cells are green, graduating through golds into salmon colors.

Map

It should be evident that an eight-point neighborhood of each cell is being used. For other neighborhoods, simply modify the nbrhood value near the beginning of expand: it is a list of index offsets relative to any given cell. For instance, a "D4" neighborhood could be specified as matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

It is also apparent that this method of spreading has its problems: it leaves holes behind. If that is not what was intended, there are various ways to fix this problem. For instance, keep the available cells in a queue so that the earliest cells found are also the earliest ones filled up. Some randomization can still be applied, but the available cells will no longer be chosen with uniform (equal) probabilities. Another, more complicated, way would be to select available cells with probabilities that depend on how many filled neighbors they have. Once a cell becomes surrounded, you can make its chance of selection so high that few holes would be left unfilled.

I will finish by commenting that this is not quite a cellular automaton (CA), which would not proceed cell by cell, but instead would update whole swathes of cells in each generation. The difference is subtle: with the CA, the selection probabilities for cells would not be uniform.

#
# Expand a patch randomly within indicator array `x` (1=unoccupied) by
# `n.size` cells beginning at index `start`.
#
expand <- function(x, n.size, start) {
  if (x[start] != 1) stop("Attempting to begin on an unoccupied cell")
  n.rows <- dim(x)[1]
  n.cols <- dim(x)[2]
  nbrhood <- matrix(c(-1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1), nrow=2)
  #
  # Adjoin one more random cell and update `state`, which records
  # (1) the immediately available cells and (2) already occupied cells.
  #
  grow <- function(state) {
    #
    # Find all available neighbors that lie within the extent of `x` and
    # are unoccupied.
    #
    neighbors <- function(i) {
      n <- c((i-1)%%n.rows+1, floor((i-1)/n.rows+1)) + nbrhood
      n <- n[, n[1,] >= 1 & n[2,] >= 1 & n[1,] <= n.rows & n[2,] <= n.cols, 
             drop=FALSE]             # Remain inside the extent of `x`.
      n <- n[1,] + (n[2,]-1)*n.rows  # Convert to *vector* indexes into `x`.
      n <- n[x[n]==1]                # Stick to valid cells in `x`.
      n <- setdiff(n, state$occupied)# Remove any occupied cells.
      return (n)
    }
    #
    # Select one available cell uniformly at random.
    # Return an updated state.
    #
    j <- ceiling(runif(1) * length(state$available))
    i <- state$available[j]
    return(list(index=i,
                available = union(state$available[-j], neighbors(i)),
                occupied = c(state$occupied, i)))
  }
  #
  # Initialize the state.
  # (If `start` is missing, choose a value at random.)
  #
  if(missing(start)) {
    indexes <- 1:(n.rows * n.cols)
    indexes <- indexes[x[indexes]==1]
    start <- sample(indexes, 1)
  }
  if(length(start)==2) start <- start[1] + (start[2]-1)*n.rows
  state <- list(available=start, occupied=c())
  #
  # Grow for as long as possible and as long as needed.
  #
  i <- 1
  indices <- c(NA, n.size)
  while(length(state$available) > 0 && i <= n.size) {
    state <- grow(state)
    indices[i] <- state$index
    i <- i+1
  }
  #
  # Return a grid of generation numbers from 1, 2, ... through n.size.
  #
  indices <- indices[!is.na(indices)]
  y <- matrix(NA, n.rows, n.cols)
  y[indices] <- 1:length(indices)
  return(y)
}
#
# Create an interesting grid `x`.
#
n.rows <- 3000
n.cols <- 5000
x <- matrix(1, n.rows, n.cols)
ij <- sapply(1:n.cols, function(i) 
      c(ceiling(n.rows * 0.5 * (1 + exp(-0.5*i/n.cols) * sin(8*i/n.cols))), i))
x[t(ij)] <- 0; x[t(ij - c(1,0))] <- 0; x[t(ij + c(1,0))] <- 0
#
# Expand around a specified location in a random but reproducible way.
#
set.seed(17)
system.time(y <- expand(x, 250, matrix(c(5, 21), 1)))
#
# Plot `y` over `x`.
#
library(raster)
plot(raster(x[n.rows:1,], xmx=n.cols, ymx=n.rows), col=c("#2020a0", "#f0f0f0"))
plot(raster(y[n.rows:1,] , xmx=n.cols, ymx=n.rows), 
     col=terrain.colors(255), alpha=.8, add=TRUE)

With slight modifications we may loop over expand to create multiple clusters. It is advisable to differentiate the clusters by an identifier, which here will run 2, 3, ..., etc.

First, change expand to return (a) NA at the first line if there is an error and (b) the values in indices rather than a matrix y. (Don't waste time creating a new matrix y with each call.) With this change made, the looping is easy: choose a random start, try to expand around it, accumulate the cluster indexes in indices if successful, and repeat until done. A key part of the loop is to limit the number of iterations in case that many contiguous clusters cannot be found: this is done with count.max.

Here is an example where 60 cluster centers are chosen uniformly at random.

size.clusters <- 250
n.clusters <- 60
count.max <- 200
set.seed(17)
system.time({
  n <- n.rows * n.cols
  cells.left <- 1:n
  cells.left[x!=1] <- -1 # Indicates occupancy of cells
  i <- 0
  indices <- c()
  ids <- c()
  while(i < n.clusters && length(cells.left) >= size.clusters && count.max > 0) {
    count.max <- count.max-1
    xy <- sample(cells.left[cells.left > 0], 1)
    cluster <- expand(x, size.clusters, xy)
    if (!is.na(cluster[1]) && length(cluster)==size.clusters) {
      i <- i+1
      ids <- c(ids, rep(i, size.clusters))
      indices <- c(indices, cluster)
      cells.left[indices] <- -1
    }                
  }
  y <- matrix(NA, n.rows, n.cols)
  y[indices] <- ids
})
cat(paste(i, "cluster(s) created.", sep=" "))

Here's the result when applied to a 310 by 500 grid (made sufficiently small and coarse for the clusters to be apparent). It takes two seconds to execute; on a 3100 by 5000 grid (100 times larger) it takes longer (24 seconds) but the timing is scaling reasonably well. (On other platforms, such as C++, the timing should hardly depend on the grid size.)

60 clusters

R – Removing Land Area from Marine Region Polygons

This might be a way.

I use this World map shapefile

# Read in new shapefile
world <- readOGR("~/test", layer="TM_WORLD_BORDERS-0.3")
# Order after $REGION in the 8th column
world_attri <- world@data[ order(-world@data[,8], world@data[,1]), ]
# Extract first line of each unique $REGION
world_attri <- world_attri[!duplicated(world_attri$REGION),]
# Subset the data to make the number of field equal to the Bahamian_shp
world_attri_subset <- subset(world_attri, select=c("NAME", "ISO3", "REGION"))
# Rename those fields to the same as the Bahamian_shp
colnames(world_attri_subset) <- c('name', 'id', 'mrgid')
# Reduce it by making a union based on $REGION
world_merged <- unionSpatialPolygons(world, world$REGION)
# Set the rownames to be equal to the mrgid field
rownames(world_attri_subset) <- world_attri_subset$mrgid
# Make a new SPDF
world_spdf <- SpatialPolygonsDataFrame(Sr=world_merged, data=world_attri_subset)

#### Following the above procedure ####

# Polygons seem to be plotted in alphabetical, so in order so make sure the bahamian ecoregion is plotted last we rename it
bahamian@data$name <- 'zBahamian'

# Insert the bahamian SPDF into the one for the world_spdf
new_spdf <- spRbind(bahamian_shp, world_spdf)
# ggplot requires a data frame, so use the fortify() function
mydf <- fortify(new_spdf, region = "name")

# Make a distinction between the underlying land shapes and the overlying ecoregion  polygon so that we can manually set the colours
mydf$filltype <- ifelse(mydf$id == 'zBahamian', "colour1", "colour2")

# Make the plot
ggplot(mydf, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", size = 1, aes(fill = mydf$filltype)) +
  + coord_equal() +
  scale_fill_manual("Test", values = c(alpha("cyan", 0.4), "white"), labels = c("Bahamian", "World"))

### Calculate the corrected area (i.e. without land area) ###

diff <- gDifference(bahamian_shp, world_spdf, byid=FALSE, id=NULL)
gArea(SpatialPolygons(diff@polygons))
# 83.69948    ## Seems legit if one compares with the below area of the whole Bahamian

gArea(SpatialPolygons(bahamian_shp@polygons))
# 84.87498

(Area in square degrees since it's not projected in UTM)

enter image description here

This approach is probably more computationally demanding compared to the previous, but at least gaining efficiency in terms of reduced manual labour.

Best Answer

Related Solutions

Raster Clumping – Creating Randomly Shaped Clumps of Cells from Seed Pixels

R – Removing Land Area from Marine Region Polygons

Related Question