Model calibration

Koen Hufkens

2026-03-31

Reading in track data

Depending on the format of your biologging data you should either read in the data from vector format (e.g. a shape file), or just a flat text file (csv, txt) with the sf package.

# reading in a text file with coordinate columns x_col and y_col
# in lon/lat format
track <- read.csv('file.csv') |>
  sf::st_as_sf(
    coords=c("x_col", "y_col"),
    crs=4326
  )

# reading in a shape file
track <- sf::st_read("your_track.shp")

Getting additional (required) data

The model framework also needs additional environmental variables to inform the behaviour of your subjects. You can use the hr_raster_maps() function to download this data. Below the routine downloads a digital elevation map (DEM) and a land cover map from the Microsoft Planetary Computer data services.

Important to note is that all settings are correctly completed or the download might fail. Consult the data services information pages to comfirm the asset name and production year (or desired year if there is a range) to download.

settings <- data.frame(
  collection = c("esa-worldcover", "cop-dem-glo-30"),
  year = c(2021, 2021),
  asset = c("map", "data"),
  filename = c("LC.tif", "DEM.tif"),
  path = "."
)

homeranger::hr_raster_maps(
  track,
  settings = settings,
  buffer = 20
)

Note that this data often needs to be converted, depending on the assumptions you make in your modelling framework. In this example the setup of Ranc et al. 2022 is used. Some conversions need to be made to get the data in to this shape.

First of all we need to read in the files, downloaded above. Both files have a different pixel resolution so we resample the lower resolution DEM data (30m) to the highest resolution (10m) the land cover product (LC). Doing so makes both raster maps compatible.

lc <- rast("LC.tif")
dem <- rast("DEM.tif")
dem <- resample(dem, lc, method = "bilinear")

We can then calculate the slope from the DEM (and its square), as well as the tree cover with a resolution of ~330m (33 pixels) using a moving window analysis. agriculture and reforested shrubland layers are also created from the land cover map.

# calculate slope and the square of it
slope <- terra::terrain(dem, v = "slope")
slope_sq <- slope^2

# calculate tree cover (class 10)
tree_cover_330 <- terra::focal(lc == 10, w=33, fun = "mean")
tree_cover_330_sq <- tree_cover_330^2

# cropland + grassland map
agriculture <- (lc == 40 | lc == 30)
reforested <- lc == 20  # shrubland as stand in for regrowth

In a final step we need to combine the different layers. The order is not of importance when the parameter ranges in the calibration step are set wide enough. When tailoring the optimization taking into account the range of the parameter values will require you to be mindful of the raster map layer order.

raster_maps <- c(slope, slope_sq, tree_cover_330, tree_cover_330_sq, agriculture, reforested)

Model calibration

To calibrate our model we now have to specify additional settings.

settings <- list(
  metric = hr_cost,
  control = list(
    sampler = "DEzs",
    settings = list(
      burnin = 10,
      iterations = 30
    )
  ),
  par = list(
    r_l = list(lower=0.0001, upper=1, init = 0.5),
    w_l = list(lower=0.0001, upper=1, init = 0.5),
    r_d = list(lower=0.0001, upper=1, init = 0.5),
    w_d = list(lower=-1, upper=-0.0001, init = -0.5),
    r_dist = list(lower=0.0001, upper=1, init = 0.5),
    w_dist = list(lower=0.0001, upper=1, init = 0.5),
    step_length_dist = list(lower=0.0001, upper=0.1, init = 0.5),
    step_length_shape = list(lower=0.3, upper=3, init = 1),
    threshold_approx_kernel = list(lower=3000, upper=10000, init = 7000),
    threshold_memory_kernel = list(lower=3000, upper=10000, init = 1000),

    # resource selection coefficients come last
    # these are unnamed
    coef = list(
      lower = rep(-3, 6),
      upper = rep(3, 6),
      init = rep(0, 6)
    )
  )
)
# calibrate the model and optimize free parameters
pars <- homeranger::hr_fit(
  drivers = raster_maps,
  obs = track,
  settings = settings,
  parallel = FALSE
)