How to create a grid to do kriging interpolation that correspond to areas in a map having actual observations?

I need to do kriging interpolation but wanted to have the interpolation in areas where I have actual data; not across the whole map.

This is my example.

Here is the base map (shp file in a .zip file): http://www.filedropper.com/lpr000b16ae

Here is the data I want to interpolate: http://www.filedropper.com/yielddata

library(dplyr)
library(ggplot2)
library(sf)
library(gstat)
library(stars)

## a) Defining the projection of the maps
proj_lambert <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

## b) Map Canada selected provinces.
canada.map <- st_read("C:/Users/rotundjo/Desktop/lpr_000b16a_e/lpr_000b16a_e.shp") %>% 
  filter(PRNAME == "Manitoba" |
           PRNAME == "Alberta" |
           PRNAME == "Saskatchewan") %>%
  st_transform(proj_lambert )

## c) Data set
data.yield <- read_csv("C:/Users/rotundjo/Desktop/yield.data.csv")
data.yield.coords <- data.yield %>% select(lat, long)
data.yield.data <- data.yield %>% select(av.yield)

#### Creating LAT LONG SpatialPoints
coordinates(data.yield.coords) = c("long", "lat")
proj4string(data.yield.coords) <- CRS("+init=epsg:4326") # 4326 is the EPSG identifier of WGS84.

### Converting the LAT LONG to the lambert CRS shown above
lambert_data.yield.coords = spTransform(data.yield.coords, proj_lambert)
lambert_data.yield.coords = as.data.frame(lambert_data.yield.coords)

#### Adding the Type column back to the data frame, with the new polar coordinates
lambert_data = cbind(lambert_data.yield.coords, data.yield.data)


coordinates(data.yield) = c("long", "lat")
proj4string(data.yield) <- CRS("+init=epsg:4326")
data.yield = spTransform(data.yield, proj_lambert)
data.yield.sf <- st_as_sf(data.yield, coords = c("long", "lat"), crs = proj_lambert)


## d) Variogram fit
v = variogram(av.yield~1, data.yield.sf)
plot(v, plot.numbers = TRUE)

v.m = fit.variogram(v, vgm("Sph", psill = 3, range = 3000000, nugget = 1))
plot(v, v.m, plot.numbers = FALSE)

## e) Building a empty grid
st_bbox(canada.map) %>%
  st_as_stars(dx = 1000) %>%
  st_set_crs(proj_lambert) %>%
  st_crop(canada.map) -> grd

## f) Kriging interpolation
k = krige(av.yield~1, data.yield.sf, grd, v.m)

## g) plotting all together
library(viridis)
ggplot() + 
  geom_stars(data = k, aes(fill = var1.pred, x = x, y = y)) +
  geom_sf(data = canada.map, color = "gray", fill="transparent") +
  geom_point(data = lambert_data, aes(x = long, y = lat), size = 0.1, color = "gray") +
  theme_bw()

library(RColorBrewer)

ggplot() + 
  geom_stars(data = k, 
             aes(fill = cut(var1.pred,breaks=c(0,20,21,22,23,24,Inf),
                            labels = c("Below 20",
                                       ">20-21", 
                                       ">21-22", 
                                       ">22-23",
                                       ">23-24",
                                       "Above 24")),
                            x = x, y = y)) +
  scale_fill_brewer(palette = "Spectral") +
  geom_sf(data = canada.map, color = "gray", fill="transparent") +
  geom_point(data = lambert_data, aes(x = long, y = lat), size = 0.1, color = "gray") +
  theme_bw()

This is my result:

But I would like to have the interpolation made only in the areas having actual observations. That is, the grd that is defined above should be defined by an area defined by the actual points. The expected result should be something like this:

The white stripped are should be part of the grd for the interpolation.

I was not able to define the grid area for interpolation to be restricted somehow to the areas where there are actual observations.