If the data is in {x,y,z}
form you can use ListDensityPlot
.
f[s_, k_] = Sin[s 2 Pi] Sin[k 2 Pi];
data = Flatten[Table[{s, k, f[s, k]}, {s, 0, 1, 0.05}, {k, 0, 1, 0.05}], 1];
ListDensityPlot[data, PlotLegends -> True]

If it is just a 2D array of data you can use MatrixPlot
f[s_, k_] = Sin[s 2 Pi] Sin[k 2 Pi];
data = Table[f[s, k], {s, 0, 1, 0.05}, {k, 0, 1, 0.05}];
MatrixPlot[data, PlotLegends -> True]

Or with discrete points
f[s_, k_] = Sin[s 2 Pi] Sin[k 2 Pi];
data = Flatten[Table[{s, k, f[s, k]}, {s, 0, 1, 0.1}, {k, 0, 1, 0.1}], 1];
min = Min@data[[All, 3]]; max = Max@data[[All, 3]];
Grid[{{Graphics[{ColorData["Rainbow"][(#[[3]] - min)/(max - min)],
PointSize[Large], Point[#[[1 ;; 2]]]} & /@ data,
ImageSize -> 300, Frame -> True], BarLegend[{"Rainbow", {min, max}}]}}]

You can choose any colorscheme you want.
RegionPlot[1 < f[s, k] < 6, {s, -2, 2}, {k, -2, 2}]
? $\endgroup$Listplot
.. because you know the function will become as in the figures in the referred post , separate points .. $\endgroup$Piecewise
.. it's only conditional function being plotted in its 2 variables plan .. $\endgroup$