Calculate the offsets (surface areas)

This is needed to fit intensities instead of numbers. The surface areas serve as weights for the mesh vertices.

dmesh_sf <- st_as_sf(dmesh)
st_crs(dmesh_sf) <- st_crs(map)
dmesh_sf <- dmesh_sf |> 
  mutate(id = row_number()) |> 
  st_sf(agr = "identity")

w <- st_intersection(dmesh_sf, st_geometry(map)) |>
  (\(x) mutate(x, area = st_area(x) |> as.numeric()))() |>
  st_drop_geometry() |>
  left_join(x = dmesh_sf, y = _, join_by(id)) |>
  mutate(area = ifelse(is.na(area), 0, area)) |>
  pull(area)