Raster R – Calculating Mean Values of All Cells in Raster Files using R

averagemeanrrasterterra

I have monthly raster (.tiff) files such as 200501.tif, 200502.tif and so on…till 202012.tif.

I want to calculate mean values of each raster (i.e. of all the cells (all rows and columns) of each raster file) (cellStats(r, 'mean')) and export them in a .csv or .xls as two columns:

Date     Value
Jan 2005  25
Feb 2005  36
...
...
...
Dec 2020  12

I am doing it in R.

Best Answer

You can use global function from {terra} library:

library(terra)

r = rast()

values(r) = 1:ncell(r)

global(r, 'mean')
#          mean
# lyr.1 32400.5

For a reproducible example I'm adding a list of random rasters loaded into a list:

library(terra)

set.seed(123)

r = rast()
rl <- sapply(1:48, function(X) setValues(r, rnorm(ncell(r))))

For computing mean for every element of the raster list, use lapply with global function:

means = lapply(rl, function(x) global(x, 'mean'))

means = do.call(rbind.data.frame, means)
row.names(means) = paste0(rep(2005:2008,each=12),sprintf("%02d",rep(1:12,times=4)))

head(means)
#                 mean
# 200501  0.0001132137
# 200502  0.0046515731
# 200503  0.0041824118
# 200504  0.0010558954
# 200505 -0.0035873522
# 200506  0.0027042357

For saving the result to your machine, use write.csv:

write.csv(means, 'path/to/means.csv')

Adaptation to OP code:

library(terra)

f = list.files(path="D:/Test/", pattern=".tiff$", full.names=TRUE)
rl = lapply(f, function(x) rast(x))

means = lapply(rl, function(x) global(x, 'mean'))
means = do.call(rbind.data.frame, means)

write.csv(means, 'path/to/means.csv')

Related Solutions

[GIS] How to perform loop on raster files

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)

R – Calculating Proportion of Polygon Covered by Each Raster Category Using R

Thank you for sending me in the right direction with exactextractr, @Spacedman. Here's the solution that I came up with:

library(terra)
library(exactextractr)
library(dplyr)
library(sf)


f <- system.file("ex/lux.shp", package="terra")
v <- st_read(f)
v <- v[1:2,]
z <- rast(vect(v), resolution=.09, names="test")
values(z) <- 1:ncell(z)
levels(z) <- 0:ncell(z)

plot(z)
plot(v, add = T)

#extract the area of each cell that is contained within each polygon
x <- exact_extract(z, v, coverage_area = TRUE)

#add polygon names that the results will be grouped by
names(x) <- v$NAME_2

#bind the list output into a df and calculate the proportion cover for each category
test <- bind_rows(x, .id = "region") %>%
  group_by(region, value) %>%
  summarize(total_area = sum(coverage_area)) %>%
  group_by(region) %>%
  mutate(proportion = total_area/sum(total_area))

Edit: With a big dataset, I was having memory issues with the code above. Here's a more memory efficient version:


### Summarizing function

sum_cover <- function(x){
  list(x %>%
    group_by(value) %>%
    summarize(total_area = sum(coverage_area)) %>%
    mutate(proportion = total_area/sum(total_area)))
  
}

#extract the area of each raster cell covered by the plot and summarize
x <- exact_extract(z, v, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)


#add plot names to the elements of the output list
names(x) <- v$NAME_2

#merge the list elements into a df
test <- bind_rows(x, .id = "Plot_buffer")

Best Answer

Related Solutions

[GIS] How to perform loop on raster files

R – Calculating Proportion of Polygon Covered by Each Raster Category Using R

Related Question