[GIS] How to perform loop on raster files

looprraster

I am new to working with spatial data in R and I have run into a few issues, that I was hoping the community could help me with.
Here is the situation:
I generated a raster file of a country with a range of probabilities. Within that country, I have state-level data – # of farms.
My goal is to randomly seed a certain number of points (farms) (# of farms is stored in CSV file, “Table2”) in each state of the country, based on the probability, indicated in the raster file. As such, I clipped my raster probability file up into 16 different raster files (the number of states in this country).

I have my code written, for what I describe above. But this code needs to be performed over each of the 16 layers – the ‘ncells’ and “TotalFarm’ value needs to be updated for each file.

install.packages("raster")  
install.packages("sp")  
install.packages("rgdal")  

setwd("~XXXXXX”)  
r <- raster("Suitability_Prob_v30.tif") #Reading in 16 rasters

r_31 <- raster("Suitability_Prob_v31.tif")  
r_32 <- raster("Suitability_Prob_v32.tif")  
r_33 <- raster("Suitability_Prob_v33.tif")  
r_34 <- raster("Suitability_Prob_v34.tif")  
r_35 <- raster("Suitability_Prob_v35.tif")  
r_36 <- raster("Suitability_Prob_v36.tif")  
r_37 <- raster("Suitability_Prob_v37.tif")  
r_38 <- raster("Suitability_Prob_v38.tif")  
r_39 <- raster("Suitability_Prob_v39.tif")  
r_310 <- raster("Suitability_Prob_v310.tif")  
r_311 <- raster("Suitability_Prob_v311.tif")  
r_312 <- raster("Suitability_Prob_v31.tif")  
r_313 <- raster("Suitability_Prob_v313.tif")  
r_314 <- raster("Suitability_Prob_v314.tif")  
r_315 <- raster("Suitability_Prob_v315.tif")

TABLE2 <- read.csv(file=“XXXXX”, header=TRUE, sep=",") #Read in cvs that has  farm #'s

Here performing ONLY on the first of the 16 raster files
Set NAs to 0

r[==-3.402823e+38] <- NA  
r[is.na(r)]=0

Calculating number of cells 

ncells<-dim(r)[1]*dim(r)[2]
ncells

Calling in the # of farms for this state, stored in CSV

TotalFarm<-sum(TABLE2[,'X30'])
TotalFarm

For point selection, set number of cells (ncells) and the number of points to be selected, (specific to each state)
Weighted by the value in the cells:

ptscell = sample(ncells, TotalFarm, prob=r[], replace=TRUE)

Distribute the points throughout the grid cells
Get the cell half-width:

hs = res(r)/2

Now find the center of those points:

centers = xyFromCell(r,ptscell)

Generate random uniform points in the cell by using the center and the
half-width/height

pts = cbind(runif(nrow(centers),centers[,1]-hs[1],centers[,1]+hs[1]),
            runif(nrow(centers),centers[,2]-hs[2],centers[,2]+hs[2]))

Print and finish:

plot(r)
points(pts)
print(pts)
write.csv(pts, file = "Points_5.txt",row.names=FALSE)

My attempts to apply this code over the series of 16 states has failed, in part, because I cannot get my 16 raster files into a “raster object” that allows me to use the rasters individually with those changes. Specifically, I made a mosaic image out of the 16 rasters , however the changes that I made on the mosaic were not made on the original raster files (only on the new mosaic object), so when I need to go back and seed the files individually with the number of farms specific to each state, the modifications that I do to the mosaic, i.e. NA’s, are not there.

All of the tutorials that I have seen perform loops on rasters by building a brick or stack. But since these are not bands which have similar extents, rather, they are neighboring states, this approach is not appropriate here.

So, how do I perform a loop on rasters, that are in the “raster” class or otherwise apply the code written above, to each raster automatically?

Best Answer

I am not sure I completely understand the question, but if it is to sample points from each state raster according to the probability found within the raster, I think this solution should do:

# Bring in libraries and set random seed
library(data.table)
library(raster)
set.seed(123)

# Bring in tifs and subset to your 16 states
tifFiles <- list.files(".", pattern=".TIF$", full.names = TRUE)
tifFiles <- paste0("./Suitability_Prob_v",c(30:39,310:315),".TIF")

# How many points do we want to sample from each state?
# I am not sure what is in your csv file, so I sampled some points from a poisson
sampPoints <- data.table(State=tifFiles, N=rpois(length(tifFiles),6))

# Function to read in raster and sample points
readRasterAndSample <- function(rasFile, count) {

    # Read in raster and deal with 0s and NAs
    ras <- raster(rasFile)
    ras[ras < 0] <- 0

    # Convert to data.table
    rasDT <- as.data.frame(ras, xy=TRUE)
    setDT(rasDT)
    setnames(rasDT, c("x","y","values"))

    # Sample non-zero points according to desired count
    # Use raster probability as sample prob
    outPoints <- rasDT[values > 0]
    outPoints <- outPoints[sample(1:.N, count, prob=values, replace=TRUE)]  

    # Jitter points randomly
    outPoints[, x := x + runif(.N, -0.5*res(ras), 0.5*res(ras))] 
    outPoints[, y := y + runif(.N, -0.5*res(ras), 0.5*res(ras))] 
    outPoints[, state := basename(rasFile)]

    # Combine into list and return
    outList <- list(ras=ras, pts=outPoints)
    return(outList)
}

# Apply function to all states
rasterList <- lapply(1:nrow(sampPoints), function(x) 
    with(sampPoints[x], readRasterAndSample(rasFile=State, count=N)))

# Collect all rasters and merge
rasterAll <- do.call(merge, lapply(rasterList, "[[", "ras"))

# Collect all state points and convert to SpatialPointsDataframe
pointDat <- rbindlist(lapply(rasterList, "[[", "pts"), fill=TRUE)
coordinates(pointDat) <- ~x+y
proj4string(pointDat) <- rasterAll@crs@projargs

# Plot to check
plot(rasterAll)
plot(pointDat, add=TRUE)
Related Question