The raster package can readily provide bilinear interpolation:
r <- rasterFromXY(df)
## due diligence, just a guess here you should check
## projection(r) <- "+proj=utm +zone=32 +datum=WGS84"
## coordinates to extract
m <- matrix(c( 3454263, 5970687), ncol = 2)
extract(r, m, method = "bilinear")
## compare with
extract(r, m, method = "simple")
See ?extract - simplest usage is a query matrix of XY coordinates in
the projection used by your raster, it will helpfully transform
queries such as a "Spatial*DataFrame" if needed, as long as both
raster x and query y have sufficient projection metadata (and it's up
to you to make sure that's set right).
(Generally building a raster from "XYZ" data is sub-optimal since
there's so much redundancy in the XY coordinates, and so much room for
things to go wrong in between. But sometimes there's no better option.