Python 3d line interpolation to increase the line resolution

I'm building a Python module to model wind turbine blades, those blades are defined as 2D profiles in X, Y space along a Z axis, I have already done that. I get a profile with the same number of points. The problem is that I want to create an STL file for that blade, and my idea is to generate a surface using the profiles, using triangulation (I don't know if this is the best solution), but the airfoils profiles are too far from each other so triangulation it's not good, so I want to increase the resolution of a blade by adding points along the Z.

in this picture you can see the position of profiles:

And here I reformat the data to connect profile points in the Z direction:

I think I can use some 3D interpolation to add n_points to each line in Z direction (above picture). To do this I've tried withsplprep and splev, here is the code that i use

import matplotlib.pyplot as plt

edge = np.array(edge)
x = edge[:, 0]
y = edge[:, 1]
z = edge[:, 2]
z_points = np.linspace(z[0], z[-1], n_points)

tck, u = splprep([x, y, z], s=0)
x_knots, y_knots, z_knots = splev(u, tck)
new_points = splev(z_points, tck)

fig2 = plt.figure(2)
ax3d = fig2.add_subplot(111, projection="3d")
ax3d.plot(new_points[0], new_points[1], new_points[2], "g")
ax3d.plot(x, y, z, "bo")
fig2.show()
plt.show()

but the result is something like this for every line, whete the spline go far fom the limits:

A sample data can be found here in this link that can be imported and viewed with the following code that produces the 2 picture:

import matplotlib.pyplot as plt

fig = plt.figure()

edges = np.load("edges.npy")
ax = fig.add_subplot(111, projection="3d")
for edge in edges:
    ax.plot(edge[:, 0], edge[:, 1], edge[:, 2], "-.")
ax.set_aspect("equal")
ax.set_axis_off()
plt.show()

It looks like you were very close. Your only mistake appears to be creating a variable z_points that varies from z[0] to z[1]. When evaluating the spline result using splev, you need to pass it a parameter vector that goes from 0 to 1.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splprep, splev

plt.close("all")

edges = np.load("edges.npy")
x = edges[:, :, 0]
y = edges[:, :, 1]
z = edges[:, :, 2]

# plotting the original edges
fig, ax = plt.subplots(subplot_kw={"projection":"3d"})
for _x, _y, _z in zip(x, y, z):
    ax.plot(_x, _y, _z, "-.")
ax.set_aspect("equal")
ax.set_axis_off()
fig.show()

# generating and plotting new edges with N points
N = 1000
t = np.linspace(0, 1, N)        # t is the parameter, which goes from 0 to 1

x_new = np.empty((x.shape[0], N))
y_new = np.empty((x.shape[0], N))
z_new = np.empty((x.shape[0], N))

fig, ax = plt.subplots(subplot_kw={"projection":"3d"})
for i, xyz in enumerate(zip(x, y, z)):
    tck, _ = splprep([*xyz], s=0)
    xyz_new = splev(t, tck)
    ax.scatter(*xyz_new, s=0.05)
    x_new[i, :], y_new[i, :], z_new[i, :] = xyz_new
ax.set_aspect("equal")
ax.set_axis_off()
fig.show()

# plotting the surface
fig, ax = plt.subplots(subplot_kw={"projection":"3d"})
ax.plot_surface(x_new, y_new, z_new, color="r", lw=0.2, edgecolor="k")
ax.set_aspect("equal")
ax.set_axis_off()
fig.show()
fig.savefig("temp.png", dpi=600, bbox_inches="tight")