I'm not familiar with this coordinate system, but I've narrowed it down to a possibility based on the info you provided. I've outlined two approaches, based on using either the "sp" or "sf" packages.
library("sp")
# As spatialpointsDataFrame
dat <- structure(list(Date = structure(c(16427, 16428, 16441, 16441,
16441, 16441, 16442, 16442, 16442, 16442, 16442, 16443, 16443,
16443, 16443, 16443, 16457, 16457, 16457, 16458), class = "Date"),
ID = c("F12", "F12", "F12", "F12", "F12", "F12", "F12", "F12",
"F12", "F12", "F12", "F12", "F12", "F12", "F12", "F12", "F12",
"F12", "F12", "F12"), X = c(238693L, 238538L, 238501L, 238507L,
238484L, 238447L, 238659L, 238507L, 238578L, 238430L, 238471L,
238534L, 238532L, 238521L, 238542L, 238580L, 238739L, 238639L,
238636L, 238759L), Y = c(2688514L, 2688467L, 2688710L, 2688718L,
2688436L, 2688421L, 2688453L, 2688702L, 2688457L, 2688452L,
2688454L, 2688440L, 2688429L, 2688698L, 2688594L, 2688621L,
2688543L, 2688507L, 2688568L, 2688497L)), class = "data.frame",
row.names = c(NA, -20L))
dat_sp <- dat
coordinates(dat_sp) <- ~X + Y
proj4string(dat_sp) <- CRS("EPSG:3826") # Details at https://epsg.io/3826
# As Simple Features object
library("sf")
dat_sf <- st_as_sf(dat, coords = c("X", "Y"), crs = 3826)
crs(dat_sf)
# To transform coordinate systems
dat_wgs84 <- st_transform(dat_sf, crs = 4326)
Coordinates plotted against an background map: