mcglinnlab / shark-ray-div

Global patterns of shark and ray diversity
0 stars 1 forks source link

Building a for loop to change resolutions #5

Closed EmmalineSheahan closed 7 years ago

EmmalineSheahan commented 7 years ago

So I'm having some trouble with this, I built this code just trying to make nine new resolutions using aggregate for factors 2 - 10 on the 110 res that's already set, but it doesn't work at all and I'm not really certain what I should actually be putting in the last space. Also this is pretty generic, you mentioned before about doing the resolution on a log base 2 scale, but I'm not really sure how to input those kinds of values. in that case I probably shouldn't use aggregate, I should just use the res() function to set a new resolution each time

factor_val <- c(2:10) new_res <- vector("list", length = length(factor_val)) for (i in seq_along(factor_val)) { new_res[[i]] <- aggregate(species_richness, fact = ) }

dmcglinn commented 7 years ago

You need to supply the argument fact with each value of factor_val you're very close! For a big block of code you can use "" in markdown (that is the text language that this issue system uses). Also you can use "r" if its R code to add some coloring to the code operations

factor_val <- c(2:10) 
new_res <- vector("list", length = length(factor_val)) 
for (i in seq_along(factor_val)) { 
    new_res[[i]] <- aggregate(species_richness, fact = factor_val[i]) 
}
dmcglinn commented 7 years ago

Oh and good call on the need to set the resolution of the raster grid rather than simply aggregate the richness raster as you have done here. The reason why this doesn't do exactly what we want is that richness is not quite additive do to species turnover. In other words if two neighboring grid cells both have 5 species then the aggregated cell may have anywhere from 5 to 10 species. See the problem?

EmmalineSheahan commented 7 years ago

Yeah, I can see that. I just made this generic for loop for changing the resolution:

factor_val <- c(200, 300, 400, 500)
new_res <- vector("list", length = length(factor_val))
new_res[[1]] <- species_richness
new_res[[2]] <- species_richness
new_res[[3]] <- species_richness
new_res[[4]] <- species_richness
for (i in seq_along(factor_val)) {
  res(new_res[[i]]) <- factor_val[i]
}

and that does work. I know there's an easier way to fill the initial values of the new_res list all with the same thing but everything I was trying was giving me an error. So I could use this to change resolutions, we just need to figure out what values we want to use for the resolutions

dmcglinn commented 7 years ago

This for loop still isn't doing quite what we want I don't think. We'll need to set the raster grid resolution for the world then rasterize each species on that grid and then sum them up. This avoids 1) polygon to grid geometry problems, 2) the non-additivity issue of richness I mentioned.

Also here is the answer to your question about initial values of a list

new_res = lapply(new_res, function(x) species_richness)

The apply functions are a family of functions for carrying out iterative tasks. In this case the function lapply is for carrying out a task on each element of a list and returning a list.

EmmalineSheahan commented 7 years ago

Ok so this changes the resolution of the world raster

factor_val <- c(110, 220, 330, 440, 550, 770, 1100)
res_list <- vector("list", length = length(factor_val))
res_list <- lapply(res_list, function(x) world_raster)
for (i in seq_along(factor_val)) {
  res(res_list[[i]]) <- factor_val[i]
}

then I'm not sure how exactly to alter the for loop for making the raster layer list, because we'd want a separate list for each different resolution

dmcglinn commented 7 years ago

Sounds like you're going to need a nested for loop! One loop to go through the different resolutions and one loop to go through the species files. For example:

for(i in 1:10) {
    for(j in 1:3) {
        print(paste("i=", i, " and j=", j, sep=''))
    }
}
EmmalineSheahan commented 7 years ago

Ok, how does this look?

# raster species polygons
sp_files = dir('./data/polygon')
sp_raster <- vector("list", length = length(sp_files))
raster_res_list <- vector("list", length = length(res_list))
for (i in seq_along(res_list)) {
     for (i in seq_along(sp_files)) {
    # read in 
    temp_poly = readOGR(dsn = paste("./data/polygon/", sp_files[i], sep=''),
                        layer = "OGRGeoJSON")
    # convert projection to cea
    temp_poly = spTransform(temp_poly, CRS("+proj=cea +units=km"))
    # add field that will provide values when rasterized
    temp_poly@data$occur = 1
    # rasterize
    sp_raster[[i]] = rasterize(temp_poly, res_list[[i]], field = 'occur')
     }
  raster_res_list[[i]] = lapply(raster_res_list, function(x) sp_raster[i])
}
dmcglinn commented 7 years ago

pretty good but the for loops cannot have the same index value in your case "i". One loop can use "i" but the other will have to use a different letter like "j" as I did in my example.

EmmalineSheahan commented 7 years ago

fixed, seems to be working thus far, just taking a little while