3d points projected to y-z plane and colorized according x-velocity using a grid for averaging

Question

This Question is related to this solved by @kglr and also this solved by @Henrik Schumacher.

I have a set of 3d positions and velocities of particles, and want to colorize the projected positions to the xz plane according to the x-velocity.

kglr solved this by usinge the following code:

SeedRandom[1];
posdata = RandomReal[10, {1000, 3}];
veldata = RandomReal[{-0.5, +0.5}, {1000, 3}];

clfun[val_] := Blend[{Blue, Green, Yellow, Red}, val];

xvelzypos = posdata;
xvelzypos[[All, 1]] = veldata[[All, 1]];

{xmin, xmax} = MinMax@veldata[[All, 1]];

legend = BarLegend[{clfun[Rescale[#, {xmin, xmax}]] &, {xmin, xmax}}, 
   LegendLabel -> "vx", LegendMarkerSize -> {20, 300}, 
   LabelStyle -> {FontFamily -> "Calibri", FontSize -> 15}];

xvzyplot = 
  Show[Graphics[{PointSize[0.03], 
       clfun[Rescale[#[[1]], {xmin, xmax}]], 
       Point[{#[[3]], #[[2]]}]} & /@ xvelzypos, Frame -> True, 
    FrameLabel -> {{"y", ""}, {"z", ""}}, AspectRatio -> Automatic, 
    PlotRange -> All, 
    BaseStyle -> {FontWeight -> "Bold", FontSize -> 20}, 
    ImageSize -> 400]];

Legended[Graphics[
  Inset[xvzyplot, Scaled[{.5, .5}], Automatic, Scaled[1]], 
  AspectRatio -> ImageAspectRatio@xvzyplot], legend]

The result is:

My question:

I would like to use a certain grid on the yz plane with values dy and dz. Each grid element should be color coded with the mean velocity of all particles that are projected to the corresponding grid element.

I tried to use parts of the solution of @Henrik Schumacher, but did not succeed.

How can be my problem solved?

Answer 1

Here is one a out of probably a million possibilities to achieve that:

SeedRandom[1];
posdata = RandomReal[10, {1000, 3}];
veldata = RandomReal[{-0.5, +0.5}, {1000, 3}];
ymin = 0.;
ymax = 10.;
zmin = 0.;
zmax = 10.;
dy = .5;
dz = .25;
m = Quotient[ymax - ymin, dy];
n = Quotient[zmax - zmin, dz];
boxcenters = Tuples[{
 Range[ymin + 0.5 dy, ymax - 0.5 dy, dy], 
 Range[zmin + 0.5 dz, zmax - 0.5 dz, dz]
 }];

If the aspect ration of the boxes is not equal to 1, the trick is to rescale the positions so that we have quadratic boxes. This way, we can employ Nearest with DistanceFunction -> ChessboardDistance to find the correct box for each point.

normalizedboxcenters = Tuples[{Range[.5, m - 0.5], Range[.5, n - 0.5]}];
idx = Flatten@Nearest[
    normalizedboxcenters -> Automatic,
    Transpose[{(#2 - ymin)/dy, (#3 - zmin)/dz} & @@ Transpose[posdata]]
    ,
    {1, 0.5},
    DistanceFunction -> ChessboardDistance
    ];

Have a index of the box for each point, we build a sparse matrix that will do the summing up by matrix-vector multiplication. Total[A, {2}] tells us, how many points are in each box which also helps us to perform the averaging.

A = SparseArray[
   Transpose[{idx, Range[Length[idx]]}] -> 1.,
   {Length[boxcenters], Length[posdata]},
   0.
   ];
meanvelocities = (A.veldata)/Clip[Total[A, {2}], {1., \[Infinity]}];
meanxvelocities = meanvelocities[[All, 1]];

Assigning some colors and plotting the result.

cols = ColorData["DarkRainbow"] /@ Rescale[meanxvelocities];
Graphics[
 Transpose[{
   cols,
   Rectangle @@@ Transpose[{
      Transpose[Transpose[boxcenters] - 0.5 {dy, dz}],
      Transpose[Transpose[boxcenters] + 0.5 {dy, dz}]
      }]
   }]
 ]

Of course, you are free to prettify the output and to add labels and stuff. But I guess, you get the idea.

Answer 2

SparseArray and "TreatRepeatedEntries"

Related Q/As :

• Mathematica fast 2D binning algorithm
• Optimising 2D binning code
• How to create a particular 5 x 5 square colour map of the mean value of data points

vvalues = xvelzypos[[All, 1]];
epsilon = 1*^-10;
indexes = 1 + Floor[(1 - epsilon) 20 Rescale[xvelzypos[[All, {2, 3}]]]];
System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries"->(Mean[{##}]&)}];
binmeansV = SparseArray[indexes -> vvalues];
System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

MatrixPlot[binmeansV, PlotRange -> {{0, 10}, {0, 10}}, 
 BaseStyle -> Opacity[.7], ColorFunction -> clfun, 
 ImageSize -> 1 -> 40, 
 FrameTicks -> {{#, #} & @ Transpose[{#, # - .5} & @Range[.5, 10.5, 1]],
    {#, #} & @ Transpose[{#, # - .5} &@Range[.5, 10.5, 1]]}]

BinLists

Using 5 bins each of size 2 for the horizontal and vertical coordinates and a single bin for velocity:

xbinspec = ybinspec = {0, 11, 2}; 
vbinspec ={-1, 1, 2};
binmeans = Flatten[Map[Mean[Last /@ #] &, 
  BinLists[RotateLeft /@ xvelzypos, xbinspec , ybinspec , vbinspec], {-3}], {3}][[1]];
matrixplot = MatrixPlot[Riffle[#, #] &[Riffle[#, #] & /@ binmeans], 
  PlotRange -> {{0, 10}, {0, 10}}, BaseStyle -> Opacity[.7], 
  ColorFunction -> clfun, ImageSize -> 1 -> 40, DataReversed -> True,
  FrameTicks -> {{#, #} &@Transpose[{#, # - .5} &@Range[.5, 10.5, 1]], 
    {#, #} &@Transpose[{#, # - .5} &@Range[.5, 10.5, 1]]}]

Add the option Epilog -> xvzyplot[[1]] to get