NAN's in the kriging predictions in R

I am trying to perform the kriging prediction of my data.

I follow a pretty standard process.
Firstly, I prep a grid to which the kriging predictions are supposed to be inserted.
I prep data that I want to interpolate.
I make sure that CRS is the same for both of the objects. Then as a next step, the variogram is selected. Lastly, we interpolate the values using kriging.

Yet, for some reason, the returned kriging object contains all predictions that are NAN.

I tried different CRS systems and I am not able to come up with a solution to this problem.

Therefore, I would like to ask for any suggestions.

Example Here:

library(tidyverse)
library(sf)
library(raster)
library(stars)
library(gstat)
library(ggmap)

# get data. using your example, we'll take spatial data from the raster package
Regions <- getData("GADM", country = "CZ", level = 1)

Regions <- Regions[Regions$NAME_1 == 'Prague', ]

Regions <- 
  Regions %>% 
  st_as_sf()

grid_spacing <- 0.005

# Create a grid for the border --------------------------------------------

polygony <- 
  st_make_grid(Regions, 
               square = T, 
               cellsize = c(grid_spacing, grid_spacing)) %>%
  st_sf()

grid <- st_intersection(polygony, Regions)

plot(grid)

grid2 <- as_Spatial(grid$geometry)
plot(grid2)

coordinates(grid2)

spdf <- grid2
grd <- makegrid(spdf, n = 3000, pretty = F)

colnames(grd) <- c('x','y')

grd %>% plot()

grd_pts <- SpatialPoints(coords = grd, 
                         proj4string=CRS(proj4string(spdf)))

# find all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[spdf, ]
grd_pts_in %>% plot()

# transform grd_pts_in back into a data frame
gdf <- as.data.frame(coordinates(grd_pts_in)) 

gdf %>% 
  mutate(value = rnorm(nrow(gdf))) %>% 
  ggplot(aes(x, y, fill =  value)) + 
  geom_tile()

grid2 %>% plot()
gdf %>% plot()  

system <- "+proj=utm +zone=19 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
system2 <- "+proj=longlat +datum=WGS84 +no_defs "



# Data to be interpolated -------------------------------------------------

df <- 
  tibble(
long = sample(coordinates(grid2)[,1], 500),
lat = sample(coordinates(grid2)[,2], 500),
val = rnorm(500, 5000, 10)
)

df %>% 
  ggplot(aes(long, lat, color = val)) + 
  geom_point()

# Kriging Predictions Itself ----------------------------------------------

resid2 <- df %>% 
  st_as_sf(coords = c("long", "lat"), crs = system2)


cutoff = 3
lzn.vgm <- variogram(val ~ 1, data = resid2, cutoff = cutoff)
plot(lzn.vgm)

nugget <- 100 # is the value of the semivariance at zero distance.
psill <- max(lzn.vgm$gamma) - nugget # is the difference between the sill and the nugget
range <- 0.2 # distance where the model first begins to flattens out.

lzn.fit.Ste <- fit.variogram(lzn.vgm, model = vgm(
  model = "Ste",
  nugget = nugget,
  sill = psill,
  range = range))

plot(lzn.vgm, lzn.fit.Ste)


coordinates(grid2)

POLY2 <- as(grid2, "Spatial")
RES <- as(resid2, "Spatial")

crs(RES) = CRS(system2)
crs(grid2)  = CRS(system2)

st_crs(grid2) == st_crs(RES)

OK =
  krige(val ~ 1, RES, grid2,
        model = lzn.fit.Ste)

OK %>% 
  as.data.frame()