2
votes

I am trying to find out the number of species per polygon within an area (America) and store this result in a table where the number of columns is the total number of species in America, and the number of rows is the total number of polygons within that area.

For this I am trying to create first an empty matrix where the results will be stored and then create a loop through all the layers (202) in a RasterBrick file called in this example “SDM.b” (each layer = distribution species (x) ), and compute zonal statistics on the raster “regionRas” which stores polygons in which I divided my area of interest. Those polygons are called “Distritos” in the raster “regionRas” and have a unique code called “CODIGO” in the raster “regionRas” please see this pic for info about those rasters.

This is the code I have:

library(maptools)
library(sp)
library(geosphere)
library(raster)


#this is the empty matrix
SpeciesRegions <- matrix(NA, nrow=87, ncol=nlayers(SDM.b))

#here I am trying to fill that matrix 
for(i in 1:nlayers(SDM.b))
{

  SpeciesRegions[,i] <- zonal(clipped.spp.diversity.mask, regionsRas,fun='sum',digits=0, na.rm=TRUE)[,1]
}

#convert the matrix data to logical (true/false) data
SpeciesRegions <- SpeciesRegions > 0

The intended output will look something like this:

             species 1  species 2 species 3
    polygon 1   FALSE   TRUE    TRUE
    polygon 2   TRUE    FALSE   FALSE
    polygon 3   FALSE   TRUE    TRUE

but this is the error I am getting in the loop

“Error in SpeciesRegions[, i] <- zonal(clipped.spp.diversity.mask, regionsRas,  : 
  number of items to replace is not a multiple of replacement length”

I would appreciate any help...

1
Thanks for the link. I have edited my question to see if that works better. I cant add data unfortunately, all my inputs are raster files which are quite big.Tac_For
Thanks for making improvements. You can still add data though. The way to do it is to use built-in data, data files easily downloaded online from stable and reputable websites, or data that is generated within your code. I have done so in the answer below. Cheers.Hack-R
Thanks so much!!! your response is also a very clear example on how to create reproducible data.Tac_For

1 Answers

2
votes

By creating reproducible data the way I described in the comments, I was eventually able to reproduce the error:

library(maptools)
library(sp)
library(geosphere)
library(raster)

# I've created reproducible data here
r     <- raster(ncols=10, nrows=10)
r[]   <- 1:ncell(r)
SDM.b <- stack(r, r, r)

clipped.spp.diversity.mask   <- raster(ncols=10, nrows=10)
clipped.spp.diversity.mask[] <- runif(ncell(r)) * 1:ncell(r)
regionsRas   <- r
regionsRas[] <- rep(1:5, each=20)


#this is the empty matrix
SpeciesRegions <- matrix(NA, nrow=87, ncol=nlayers(SDM.b))

#here I am trying to fill that matrix 
for(i in 1:nlayers(SDM.b)){
  SpeciesRegions[,i] <- zonal(clipped.spp.diversity.mask, regionsRas,fun='sum',digits=0, na.rm=TRUE)[,1]
}
Error in SpeciesRegions[, i] <- zonal(clipped.spp.diversity.mask, regionsRas,  : 
  number of items to replace is not a multiple of replacement length

This error message can usually be taken pretty literally, so I checked the dimensions of the existing and replacement values:

zonal(clipped.spp.diversity.mask, regionsRas,fun='sum',digits=0, na.rm=TRUE)[,1]

[1] 1 2 3 4 5

SpeciesRegions[,i]
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[65] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

As you can see, they are very different lengths and thus it should be unsurprising that I got that error.

To fix it, I can change the hard-coded row length in the definition of the empty dataframe to the length of the input data:

#this is the empty matrix
SpeciesRegions <- matrix(NA, nrow=length(zonal(clipped.spp.diversity.mask, regionsRas,fun='sum',digits=0, na.rm=TRUE)[,1]), ncol=nlayers(SDM.b))

#here I am trying to fill that matrix 
for(i in 1:nlayers(SDM.b)){
  SpeciesRegions[,i] <- zonal(clipped.spp.diversity.mask, regionsRas,fun='sum',digits=0, na.rm=TRUE)[,1]
}

Now the code runs fine.