You can use scipy.interpolate
for this.
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splev, splrep
x = np.linspace(0, 10, 10) # x-coordinates
y = np.sin(x) # y-coordinates
tck = splrep(x, y) # get bspline representation given (x,y) values
x2 = np.linspace(0, 10, 200) # new set of values, just to check
y2 = splev(x2, tck) # evaluate the y values of new coordinates on NURBS curve
plt.plot(x, y, 'o', x2, y2)
plt.show()

The tuple tck
contains your knot vector and control points (coefficients). There are also more involved routines in SciPy, look here.
Note that these are only for bspline curves. As far as I know there are no equivalent methods for surfaces in SciPy. If you want to use surfaces, depending on your requirement, you can either use igakit
from igakit.cad import ruled, circle
c1 = circle(angle=(0,np.pi/2.))
c2 = circle(radius=2,angle=(0,np.pi/2.))
print "knot vector:", c1.knots
print "control points:", c1.control
srf = ruled(c1,c2)
plt.plot(srf)
plt.show()
knot vector: (array([ 0., 0., 0., 1., 1., 1.]),)
control points: array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
[ 7.07106781e-01, 7.07106781e-01, 0.00000000e+00, 7.07106781e-01],
[ 2.22044605e-16, 1.00000000e+00, 0.00000000e+00, 1.00000000e+00]])

or the NURBS package. For fancier stuff Blender and Salome have complete Python API for all family of NURBS curves/surfaces with the latter being based on OpenCascade.