nurbs_surface.py

import bpy
from mathutils import Vector

surface_data = bpy.data.curves.new('wook', 'SURFACE')
surface_data.dimensions = '3D'
spline_0 = surface_data.splines.new(type='NURBS')

# 16 coordinates
points = [
    Vector((-1.5, -1.5, 0.0, 1.0)), Vector((-1.5, -0.5, 0.0, 1.0)),
    Vector((-1.5, 0.5, 0.0, 1.0)), Vector((-1.5, 1.5, 0.0, 1.0)),
    Vector((-0.5, -1.5, 0.0, 1.0)), Vector((-0.5, -0.5, 1.0, 1.0)),
    Vector((-0.5, 0.5, 1.0, 1.0)), Vector((-0.5, 1.5, 0.0, 1.0)),
    Vector((0.5, -1.5, 0.0, 1.0)), Vector((0.5, -0.5, 1.0, 1.0)),
    Vector((0.5, 0.5, 1.0, 1.0)), Vector((0.5, 1.5, 0.0, 1.0)),
    Vector((1.5, -1.5, 0.0, 1.0)), Vector((1.5, -0.5, 0.0, 1.0)),
    Vector((1.5, 0.5, 0.0, 1.0)), Vector((1.5, 1.5, 0.0, 1.0))
]

spline_0.points.add(15)  #   already has a a default zero vector

for p, new_co in zip(spline_0.points, points):
    p.co = new_co

# 4*4 = 16 points
# spline_0.point_count_u = 4   # read only :(
# spline_0.point_count_v = 4   # read only :(
spline_0.order_v = 4
spline_0.order_u = 4
spline_0.resolution_v= 4
spline_0.resolution_u = 4

surface_object = bpy.data.objects.new('NURBS_OBJ', surface_data)

scene = bpy.context.scene    
scene.objects.link(surface_object)