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)