Spatial interpolation of climate data using latitude, longitude and elevation in r

I have downloaded E-OBS gridded daily climate data on 0.1-degree grid. For each grid point, I have longitude, latitude, elevation and climate data, i.e. precipitation sum / mean temperature. Now, I would like interpolate climate data for specific site, considering lon, lat and elevation.

So, I have achieved this by using simple (multiple) linear regression (see the code below), but by reading the relevant literature, I can inform myself that for such purposes, kriging or cokriging is usually applied, e.g. here or here.
Can someone help me with the code development, where I can apply (co)kriging on my data or maybe suggest even more sophisticated method?

 # 1 Construct a sample dataset
structure(list(lon = c(15.35, 15.45, 15.55, 15.65, 15.75, 15.35, 
                       15.45, 15.55, 15.65, 15.75, 15.35, 15.45, 15.55, 15.65, 15.75, 
                       15.35, 15.45, 15.55, 15.65, 15.75, 15.35, 15.45, 15.55, 15.65, 
                       15.75), lat = c(46.45, 46.45, 46.45, 46.45, 46.45, 46.55, 46.55, 
                                       46.55, 46.55, 46.55, 46.65, 46.65, 46.65, 46.65, 46.65, 46.75, 
                                       46.75, 46.75, 46.75, 46.75, 46.85, 46.85, 46.85, 46.85, 46.85
                       ), t_avg = c(5.39999987930059, 6.52999985404313, 9.01999979838729, 
                                    11.0299997534603, 11.169999750331, 9.18999979458749, 9.49999978765845, 
                                    9.89999977871776, 11.1599997505546, 11.0299997534603, 10.049999775365, 
                                    10.1899997722358, 10.5699997637421, 11.1599997505546, 11.1099997516721, 
                                    11.6699997391552, 11.2799997478724, 11.7199997380376, 11.6799997389317, 
                                    11.5099997427315, 12.059999730438, 11.9399997331202, 11.8699997346848, 
                                    11.499999742955, 11.2499997485429), elevation = c(1149.61, 994.63, 
                                                                                      639.31, 293.11, 245.52, 661.85, 608.48, 534.67, 296.02, 283.24, 
                                                                                      564.26, 527.31, 444.78, 318.32, 305.5, 339.65, 386.74, 298.66, 
                                                                                      271.16, 263.38, 318.88, 323.87, 293.57, 333.73, 323.76), p_sum = c(40.1000005975366, 
                                                                                                                                                         40.7000006064773, 29.0000004321337, 24.9000003710389, 24.5000003650784, 
                                                                                                                                                         35.5000005289912, 35.2000005245209, 33.8000005036592, 28.6000004261732, 
                                                                                                                                                         29.9000004455447, 26.3000003919005, 28.400000423193, 28.0000004172325, 
                                                                                                                                                         26.0000003874302, 28.1000004187226, 19.1000002846122, 20.8000003099442, 
                                                                                                                                                         20.3000003024936, 21.4000003188848, 25.1000003740191, 16.7000002488494, 
                                                                                                                                                         17.8000002652407, 19.1000002846122, 21.7000003233552, 25.1000003740191
                                                                                      )), row.names = c(NA, -25L), class = "data.frame")

# This is location for the prediction
pred_df <- data.frame(
  lat = 46.68,
  lon = 15.55,
  elevation = 210
)

# 1 Predictions for temperature
lin_model_T <- lm(t_avg ~ lon + lat + elevation, data = df_small)  
pred_df$t_avg <- predict(lin_model_T, pred_df)

# 1 Predictions for precipitation
lin_model_P <- lm(p_sum ~ lon + lat + elevation, data = df_small)
pred_df$p_sum <- predict(lin_model_P, pred_df)