1

I have a csv file of discrete point locations containing lat and long for different genera, with multiple entries for most genera, and I want to create a raster stack with a raster layer for each genus, with the rasters indicating the presence/absence of each genus. Because there are over 1700 genera in my file, I need to create some kind of loop to do this, unless there's an easier way, but for the life of me, I can't figure out how to do this. This is the code I wrote- I am aware that this is the wrong way to write loops, but I don't know what else to do. Here is my code:

cell_size <- 2.5
lon_min <- -179.0816; lon_max <- 178.5042; lat_min <- 24.5542; lat_max <- 81.817
ncols <- ((lon_max - lon_min)/cell_size)+1; nrows <- ((lat_max - lat_min)/cell_size)+1

for (i in cropaphidlocSPdata$Genus) {
  while !(i %in% aphidlist) {
    [i]data <- subset(cropaphidlocSPdata, Genus == [i]
    [i]count <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, 
                       res=2.5, crs="+proj=longlat +datum=WGS84")
    [i]coords <- cbind([i]data$Lon, [i]data$Lat)
    [i]count <- rasterize([i]coords, [i]count, fun = "count")
    aphidcounts <- addLayer(aphidcounts, [i]count)
    aphidlist <- c(aphidlist, [i])
  • Do you have multiple columns for single species/genus or are counts for a genus contained in each column (eg., only one column for a genus)? If there are multiple columns representing a given species/genus then the first step is to aggregate the data in some way. A sample of your data would be ideal in providing a reproducible example. This type of data can take many forms and you are leaving us guessing. – Jeffrey Evans Jan 25 at 20:42
  • I added to my answer to hopefully address the structure of your data. – Jeffrey Evans Jan 25 at 21:31
1

If you do not hit memory issues, the easiest way would be to use the sp raster class. If the data is, in fact, gridded then you can use read.csv to read the data then coerce into a SpatialPixelsDataFrame. This format is just like other sp class objects, containing an @data slot with the raster attributes as columns in a data.frame. To get a raster class object from this sp object, containing all of the species data, you would just use raster::stack.

library(sp)
library(raster)

data(meuse.grid)
m = SpatialPixelsDataFrame(points = meuse.grid[c("x", "y")], 
                           data = meuse.grid)
m@data <- data.frame(spp1=round(runif(nrow(m), 0, 1),0),
                     spp2=round(runif(nrow(m), 0, 1),0),
                     spp3=round(runif(nrow(m), 0, 1),0),
                     spp4=round(runif(nrow(m), 0, 1),0)) 
head(m@data)

spp <- stack(m)
  plot(spp)

Here is an example using irregular points where you have a data.frame of [x,y] locations with a single column indicating a genus at a given location. Here, I create some dummy data to approximate your data (eg., [x,y,genus]).

library(sp)
library(raster)

Create empty reference raster

r <- raster(extent(-179.0816, 178.5042, 24.5542, 81.817), resolution = 2.5) 

Create some dummy data with x,y locations and a column indicating genus presence.

genera <- data.frame(coordinates(spsample(as(extent(r), "SpatialPolygons"), 
                     n=200, type="random")), genus = sample(c("genus1", "genus2", 
                     "genus3", "genus4"), 200, replace = TRUE))                                
head(genera)

Coerce to SpatialPointsDataFrame

coordinates(genera) <- ~x+y
  plot(genera)                   

Loop to rasterize points by subset of each genus, this will return number of observations, of a genus, within each cell. If you want a binomial just change the fun argument in the rasterize function to "sum".

rs <- stack(r)                   
  for(i in unique(genera$genus)) {
      rs <- addLayer(rs, rasterize(x = genera[genera$genus == i,], y = r, 
                     field = 1, fun=function(x, ...)length(x), 
                     background = 0))
  }
names(rs) <- unique(genera$genus) 

plot(rs) 
  • What do you mean 'if the data is gridded', and why do you write it as .grid ? – birdoptera Jan 25 at 19:36
  • As far as I can tell, the data is not gridded – birdoptera Jan 25 at 19:44
0

I was able to revise my for loop into something that worked;

cell_size <- 2.5
lon_min <- -179.0816; lon_max <- 178.5042; lat_min <- 24.5542; lat_max <- 81.817
ncols <- ((lon_max - lon_min)/cell_size)+1; nrows <- ((lat_max - lat_min)/cell_size)+1


for (i in levels(cropaphidlocSPdata$Genus)) {
  print(i)
    x <- subset(cropaphidlocSPdata, Genus == i)
    ra <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, 
                       res=0.25, crs="+proj=longlat +datum=WGS84")
    coords <- cbind(x$Lon, x$Lat)
    oc <- rasterize(coords, ra, fun = "count")
    aphidcounts <- addLayer(aphidcounts,oc)

  }
  • Did my answer not work? Why are you recreating the raster in every loop, it is just a reference for rasterize? You really did not take any advice here. Nesting a good skill to learn in R, the second approach I show in my answer does all of this with one call: rs <- addLayer(rs, rasterize(x = genera[genera$genus == i,], y = r, field = 1, fun=function(x, ...)length(x), background = 0)) – Jeffrey Evans Jan 28 at 19:25
  • I was shown this approach before I read your response, so I just went with it. – birdoptera Jan 28 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.