3D binning keeping the values

Question

I have N points in the from of triplets {xi,yi,zi}.

I would like to bin them on a grid nbin^2 based on their x, y values. Now, the key difference compared to this question is that I would like to compute average values in each bin for binned number of points and obtain {xav,yav,zav} for bin 1, up to nbin^2 and to plot nbin^2 number of points instead of initial N points.

I do not want to plot it against bin values as in example above, and please note that x, y will not be equidistant. I started the code for pair {xi,yi} and it does what I want (please see below), I am stuck how to effectively extend it in 3D and ideally optimise in line with something like this.

test = RandomInteger[{-100, 100}, {1000, 2}];

xmin = Min@test[[1 ;;]][[All, 1]];
xmax = Max@test[[1 ;;]][[All, 1]];
ymin = Min@test[[1 ;;]][[All, 2]];
ymax = Max@test[[1 ;;]][[All, 2]];
nbin = 15.;

t2 = BinLists[
 test[[1 ;; ;; 1]], {xmin, xmax, (xmax - xmin)/nbin}, {ymin, 
ymax, (ymax - ymin)/nbin}];
Length /@ {t2, t2[[1]]}

(* getting rid of empty bins*)

  ct2 = Replace[t2, x_List :> DeleteCases[x, {}], {0, Infinity}];   

  ave1 = Map[Mean, ct2[[#]], 1] & /@ Range[Length[ct2]];
  ListPlot /@ {test, ave1}

Answer 1

n = 2;
nbins = 15;
testb = RandomInteger[{-100, 100}, {1000, n}];
epsilon = 1*^-10;
indexes = 1 + Floor[(1 - epsilon) nbins Transpose[Rescale /@ Transpose[testb][[{1, 2}]]]];

Make the data entries strings, and modify the "TreatRepeatedEntries" suboption value to make the list of repeated entries a string:

values = ToString /@ testb;

System`SetSystemOptions["SparseArrayOptions" ->{"TreatRepeatedEntries" ->
  (ToString @ {##} &)}];
binmeansZ = SparseArray[indexes -> values];
System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

Process the sparse array to change strings to expressions and take the means:

newlist = Replace[binmeansZ["NonzeroValues"] /. x_String :> ToExpression[x], 
   y : {ConstantArray[_, Dimensions[testb][[2]]] ..} :> Mean[y], 1];

Display:

ListPlot[{testb, newlist}, 
 PlotStyle -> {Directive[PointSize[Small], Blue], Directive[PointSize[Medium], Red]}, 
 AspectRatio -> 1, Frame -> True, Axes -> False, 
 GridLines -> {Range[-100, 100, 200/15], Range[-100, 100, 200/15]}]

For n = 3 we get

c = {Red, Blue};
ListPointPlot3D[#, ImageSize -> 300, BoxRatios -> 1, Axes -> False, 
    PlotStyle -> Directive[PointSize[Medium], First[c = RotateLeft[c]]]] & /@ 
  {testb, newlist} // Row