Automate R code

[nikos_geo]

I have 31 raster layers and I want to measure the root mean square error (RMSE) of a linear model. Each time I am using the independent variable and one of the dependent variables. The names:

independent variable: lj.tif
dependent variables: ntl_atprk000.tif, ntl_atprk010.tif, ntl_atprk020.tif...., ntl_atprk300.tif

Now, I am calculating the RMSE manually, by changing the name of the dependent variable. For example:

library(terra)

wd <- "path/"

ntl = rast(paste0(wd, "ntl_atprk000.tif")) # or ntl_atprk010.., ntl_atprk180,..., ntl_atprk300
covar = rast(paste0(wd, "lj.tif"))

covar = resample(covar, ntl, method = "bilinear")

s = c(ntl, covar)
names(s) = c("ntl", "covar")

m <- lm(ntl ~ covar, data = as.data.frame(s), na.action = na.omit)

# rmse lm
sqrt(mean(m$residuals^2))

How can I change the value of y in the linear model automatically?

For simplicity, here are three raster layer (the independent raster and two dependent rasters):

independent var:
covar = rast(ncols=447, nrows=328, nlyrs=1, xmin=-31952.684, xmax=12747.316, ymin=6012571.3953, ymax=6045371.3953, names=c('lj'), crs='PROJCRS[\"World_Mollweide\",BASEGEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"Degree\",0.0174532925199433]]],CONVERSION[\"unnamed\",METHOD[\"Mollweide\"],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"Degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]')

dependent 1 (ntl_atprk000):
ntl = rast(ncols=437, nrows=318, nlyrs=1, xmin=-31400, xmax=12300, ymin=6013100, ymax=6044900, names=c('layer'), crs='PROJCRS[\"unknown\",BASEGEOGCRS[\"GCS_unknown\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"Degree\",0.0174532925199433]]],CONVERSION[\"unnamed\",METHOD[\"Mollweide\"],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"Degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]')

dependent 2 (ntl_atprk170.tif):
ntl = rast(ncols=437, nrows=318, nlyrs=1, xmin=-31400, xmax=12300, ymin=6013100, ymax=6044900, names=c('layer'), crs='PROJCRS[\"unknown\",BASEGEOGCRS[\"GCS_unknown\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ID[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"Degree\",0.0174532925199433]]],CONVERSION[\"unnamed\",METHOD[\"Mollweide\"],PARAMETER[\"Longitude of natural origin\",0,ANGLEUNIT[\"Degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]')

[nirgrahamuk]

Step 1 is learning how to write r functions. This allows you to vary an input and get associated output

Step 2 is iterating either with r builtins like the apply family of functions or their tidyverse equivalents in purrr ; i.e. you construct the parameters against which you want your function run, and then run them together.

[nikos_geo]

The solution:

library(terra)

wd <- "path/"

covar = rast(paste0(wd, "lj.tif"))

rlist = list.files(path = wd,
                   pattern = "^ntl_atprk\\d+\\.tif$",
                   all.files = T,
                   full.names = F)

for (i in rlist){
  for (j in i) {
    
    nameNum = gsub("\\D+","",j)
    print(nameNum)
    print(j)
    
    atprk = rast(paste0(wd, j))
    
    covar = resample(covar, atprk, method = "bilinear")
    
    s = c(atprk, covar)
    names(s) = c("atprk", "covar")
    
    m <- lm(atprk ~ covar, data = as.data.frame(s), na.action = na.omit)
    
    # rmse lm
    rmse = sqrt(mean(m$residuals^2))
    
    print(rmse)
  }
}