CODEWITHSUNDEEP
rr-rasterterra
Shakoo
Shakoo

Reputation: 1

Calculate Hurst Exponent on Raster stack In R

I have a stack of five maps showing average NDVI values, and I want to calculate the Hurst exponent for each pixel in the stack using the "hurstexp" function from the "pracma" package. My stacked maps have some missing values (NAs) represented as -Inf when extracted, and I want to return an output raster with the calculated Hurst exponent values. Can you please help me with this, including how to handle the NAs in the process?


library(sp)
library(raster)
library(pracma)

modis_list <- list.files(path = "path",
                         pattern = "*.tif$", full.names = TRUE)
modis_raster <- stack(modis_list)

# i try this but not working
modis_hurst <- raster::calc(modis_raster, function(x) {
  spec <- hurstexp(x, d = 50, display = TRUE)
  return(spec$H)
})


writeRaster(modis_hurst, filename = "path/to/output/file.tif", format = "GTiff")


Error:Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

#### I also tried this
calculate_hurst <- function(x) {
  tryCatch({
    spec <- hurstexp(x, d = 50, display = FALSE)  # Use display = FALSE to suppress plot
    return(as.numeric(spec$H))  # Ensure the result is numeric
  }, error = function(e) {
    return(NA)  # Handle NA values or errors during calculation
  })
}
# Calculate the Hurst exponent for each layer using calc
modis_hurst <- raster::calc(modis_raster, calculate_hurst)

# Write the result to a new raster file
writeRaster(modis_hurst, filename = "path/to/output/file.tif", format = "GTiff")


Error in setValues(out, x) : values must be numeric, logical or factor

Upvotes: 0

Views: 204

Answers (1)

Elia
Elia

Reputation: 2584

You can use the function app in the terra package (i.e., the new raster package). Since I don't know the function hurstexp and its results meaning, I suppose that whenever a pixel in your time series is NA, the results should be NA (but this is up to you). In any case, hurstexp cannot handle NA, so you have to get rid of them.

library(terra)
library(pracma)

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
#create raster stack
rr <- c(r,r*2,r/2,r^2,log(r))
plot(rr)

spec <- terra::app(rr,function(x){
  x <- as.numeric(x)
  x <- x[!is.na(x)]
  if(length(x)==0)return(NA)
  spec <- hurstexp(x,display = F)
  Hs <- as.numeric(spec$Hs)
  return(Hs)
},cores=1)

plot(spec)

If you (likely) have more than 1 core in your machine you can change the cores argument of app to speed up the computation by running the function in parallel.

EDIT

The problem with returning more than one parameter in hutstexp is that app cannot return a raster stack. So, one solution could be to loop over the desired parameters with e.g. lapply.

res <- lapply(list("Hs","He"),function(par){
  spec <- terra::app(rr,function(x){
    x <- as.numeric(x)
    x <- x[!is.na(x)]
    if(length(x)==0)return(NA)
    spec <- hurstexp(x,display = F)
    res <- spec[[par]]
    return(res)
  },cores=1)
})

res <- do.call(c,res)
plot(res)

Maybe there are better ways, but an advantage of this approach is that you can parallelize both the function in app with the argument cores and the outer lapply loop, with e.g., foreach or snowfall (but all other parallelization mechanisms should work fine)

Upvotes: 0

Related Questions