Overlapping splot colourmap with 2-d points Gnuplot

It seams a simple question but I didn't find a solution elsewhere. I am plotting a colourmap with splot and would like to add to the plot the points at which the space was sampled.

I'll show the code using a 3 dimensional case, but it could be 3-d or higher

Set and unset used:

number_of_variables = 3
set contour base
set palette model RGB
set palette defined ( 0 'white', 1 'dark-red' )
set dgrid3d 31, 31
unset clabel
unset surface
set view map
set pm3d 
set datafile missing NaN

Section where the splot is used:

do for [i=1:number_of_variables-1]{ 
  set xlabel sprintf("b %d", i)
  do for [j=i+1:number_of_variables] {
    set ylabel sprintf("b %d", j)
    set title sprintf("Kriging expected improvement with sampling point %d vs %d", i, j)
    splot "sf1.dat" using i:j:(((column(number_of_variables+1) == (i-1))||(column(number_of_variables+1) == -1)&&((column(number_of_variables+2) == (j-1))||(column(number_of_variables+2) == -1))) ? column(number_of_variables+5) : 1/0) notitle w l lt -1, \
            "sf2.dat" using i:j:(0) w p pt 7 ps 1.5 lc rgb "black" title "Sampled points"
    pause -1 
  }
}

Example data for a 3d case:

"$sf1.dat" has 6 columns: x y z of each point, plotted 2 dimensions at a time, 2 columns used in the rule and 1 column for the value of the colourmap number of rows are 31^3

-0.150764 -0.0776023 0.00485263 0 1 0.333911
-0.150764 0.066001 0.0012847 0 1 0.343821
etc...

"$sf2.dat" has 3 columns: x y z and I want to plot 2 dimensons at a time number of rows are around 35

0.77489 0.866638 0.697351 
0.175989 0.00346154 1.01085 
etc... 

The rule for the third column chooses if to pass the last column or 1/0 based on the value of the 4th and 5th column for the case above.

The problem:

splot "sf1.dat" using i:j:rule notitle w l lt -1,\
      "sf2.dat" using i:j:(0) w p pt 7 ps 1.5 lc rgb "black" title "Sampled points"

Where splot is used, I get a white rectangle on the graph (white because I assign value 0) instead of the dots at the sampled locations.

The output I get - image

Image of sf1.dat I want to plot the dots on

If I use a table and plot both in 2D I get the dots at the correct location, and instead of the colourmap it shows the contour lines (I would like the opposite).

Full data: sf1.dat sf2.dat

If I missed anything please add a comment. Thank you for the help

Edit 1

I am using Gnuplot verion 5.4 patchlevel 2, last modified 01-06-2021.

I'm not 100% sure whether the following is the way you want it to look like. My first thought would be to plot the data via dgrid3d into a table and then plot it as 2D-plot and overlaying the other points. This should work for gnuplot>=5.0. I read somehwere that gnuplot>=5.4.3 might have some options to do it without plotting to a table first, but since you have 5.4.2...
You can also play with qnorm <value>, check help dgrid3d.

Script: (works for gnuplot>=5.4.0, apparently needs some adaptions for <5.4.0)

### point data on top of dgrid3d data
reset session

FILE1 = "SO79089934_sf1.dat"
FILE2 = "SO79089934_sf2.dat"

N       = 3    # number_of_variables
rule(n) = ((column(n+1) == (i-1))||(n+1) == -1) && \
          ((column(n+2) == (j-1))||(column(n+2) == -1)) ? column(n+5) : 1/0

set palette defined ( 0 'white', 1 'dark-red' )
set size ratio -1
set xrange noextend
set yrange noextend
set key noautotitle
unset contour
set dgrid3d 31,31 qnorm 1

do for [i=1:N-1]{ 
    do for [j=i+1:N] {
        set xlabel sprintf("b %d", i)
        set ylabel sprintf("b %d", j)
        set table $GridFile1
            splot FILE1 u 1:2:(rule(N))
        unset table
        set title sprintf("Kriging expected improvement with sampling point %d vs %d", i, j)

        plot $GridFile1 u 1:2:3 w image lc palette z, \
                  FILE2 u i:j w p pt 7 ps 1.5 lc rgb "black" ti "Sampled points"
        if (i*j < N*(N-1)) { pause -1 }
    }
}
### end of script

Result:

This one gives a warning:
Warning: empty x range [1.42008:1.42008], adjusting to [1.40588:1.43428]