Color of a point-made 3D surface

I have a set of 4D data (3D points + 1D value). Let's call them X,Y,Z and C. I would like to generate a surface from the X,Y,Z points and then colour it according to the C values.<br> I think I'm asking the same thing that <i>diffracteD</i> did in this question, but nobody seems to have understood what he was asking and the answers and comments aren't helpful.

I was able to create a surface with XYZ data, following this this answer, but then the surface is coloured according to the Z value, not the C value as I want. <br>

On the other end, the <i>original answer</i> to this question manages to colour the surface using the C value, but in this case Z is a function of X and Y, not a free variable as in my case.

My goal is to somewhat merge the two things, creating a surface from XYZ data and colouring it according to C

Here is what I've done so far:

1. XYZ independent but surface coloured with Z values:

import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import cm

import numpy as np
from numpy.random import randn
from scipy import array, newaxis


# ======
## data:

DATA = np.array([
    [-0.807237702464, 0.904373229492, 111.428744443],
    [-0.802470821517, 0.832159465335, 98.572957317],
    [-0.801052795982, 0.744231916692, 86.485869328],
    [-0.802505546206, 0.642324228721, 75.279804677],
    [-0.804158144115, 0.52882485495, 65.112895758],
    [-0.806418040943, 0.405733109371, 56.1627277595],
    [-0.808515314192, 0.275100227689, 48.508994388],
    [-0.809879521648, 0.139140394575, 42.1027499025],
    [-0.810645106092, -7.48279012695e-06, 36.8668106345],
    [-0.810676720161, -0.139773175337, 32.714580273],
    [-0.811308686707, -0.277276065449, 29.5977405865],
    [-0.812331692291, -0.40975978382, 27.6210856615],
    [-0.816075037319, -0.535615685086, 27.2420699235],
    [-0.823691366944, -0.654350489595, 29.1823292975],
    [-0.836688691603, -0.765630198427, 34.2275056775],
    [-0.854984518665, -0.86845932028, 43.029581434],
    [-0.879261949054, -0.961799684483, 55.9594146815],
    [-0.740499820944, 0.901631050387, 97.0261463995],
    [-0.735011699497, 0.82881933383, 84.971061395],
    [-0.733021568161, 0.740454485354, 73.733621269],
    [-0.732821755233, 0.638770044767, 63.3815970475],
    [-0.733876941678, 0.525818698874, 54.0655910105],
    [-0.735055978521, 0.403303715698, 45.90859502],
    [-0.736448900325, 0.273425879041, 38.935709456],
    [-0.737556181137, 0.13826504904, 33.096106049],
    [-0.738278724065, -9.73058423274e-06, 28.359664343],
    [-0.738507612286, -0.138781586244, 24.627237837],
    [-0.738539663773, -0.275090412979, 21.857410904],
    [-0.739099040189, -0.406068448513, 20.1110519655],
    [-0.741152200369, -0.529726022182, 19.7019157715],
])

Xs = DATA[:,0]
Ys = DATA[:,1]
Zs = DATA[:,2]

# ======
## plot:

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')


#cmap="hot" colours w.r.t. Z values

surf = ax.plot_trisurf(Xs, Ys, Zs, cmap="hot", linewidth=0)
fig.colorbar(surf)

ax.xaxis.set_major_locator(MaxNLocator(5))
ax.yaxis.set_major_locator(MaxNLocator(6))
ax.zaxis.set_major_locator(MaxNLocator(5))

fig.tight_layout()

plt.show()

2. Surface coloured with C values, but Z is a function of X,Y

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

# as plot_surface needs 2D arrays as input
x = np.arange(10)
y = np.array(range(10,15))
# we make a meshgrid from the x,y data
X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X**2 + Y**2))

# data_value shall be represented by color
data_value = np.random.rand(len(y), len(x))
# map the data to rgba values from a colormap
colors = cm.ScalarMappable(cmap = "hot").to_rgba(data_value)


# plot_surface with points X,Y,Z and data_value as colors
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=colors,
                       linewidth=0, antialiased=True)
fig.colorbar(surf)
plt.show()

Is there a way to modify the first code to make it use C values instead of Z values to colour the surface?
Thanks.

Here is an alternative method with setting the face colors to the triangles in the trisurf, just need to make the C values have the same number of the triangles:

Xs = DATA[:,0]
Ys = DATA[:,1]
Zs = DATA[:,2]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_trisurf(Xs, Ys, Zs, linewidth=0)

We can find the triangles in the figure clearly :

Then just need to generate the colors for the triangles:

# Find the slices of the triangles
slices = surf._segslices

# data_value shall be represented by color
data_value = np.random.rand(len(slices))
# map the data to rgba values from a colormap
colors = cm.ScalarMappable(cmap = "hot").to_rgba(data_value)

# set the face colors
surf.set_fc(colors)
cbar = fig.colorbar(cm.ScalarMappable(cmap="hot"), ax=ax)

ax.xaxis.set_major_locator(MaxNLocator(5))
ax.yaxis.set_major_locator(MaxNLocator(6))
ax.zaxis.set_major_locator(MaxNLocator(5))
fig.tight_layout()

plt.show()

Get the figure :

I ended up following <i>HMH1013</i> suggestion and use plotly.
Here is the code:

import plotly.graph_objects as go
import numpy as np

#data.dat has 4 cols of numbers: X, Y, Z, C
DATA = np.loadtxt(open("data.dat", "rb"))

Xs = DATA[:,0]
Ys = DATA[:,1]
Zs = DATA[:,2]
values = DATA[:,3]

fig = go.Figure(data=[
    go.Mesh3d(
        x=Xs,
        y=Ys,
        z=Zs,
        colorbar_title='value',
        colorscale="jet",
        intensity=values,
        showscale=True
    )
])

fig.show()

In this way C needs to have the same dimension of X,Y and Z and not the one of the triangles.

Many thanks to <i>HMH1013</i> for helping me.