1

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?

  • Take a look at stack and stackApply. You can index a raster stack using a double bracket, eg., r[[i]] making looping over a stack quite easy. – Jeffrey Evans Jan 24 '17 at 2:51
1

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)

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.