Need to fit a 3D surface over a grid of 3D data

[Roman Rodriguez]

My intention is to fit 3-D data (x, y, z), where z has no known function f(x,y) -> z. To be precise, I want to create a volatility surface from an option chain for S&P500, where the data is (x, y, z) -> (time, option strikes, volatility). While it is true that implied volatility depends somehow on time and strike, it really depends on many more parameters, so let's assume the f(x,y) function is unknown.

Therefore, if I just have a bunch of 3D data, I need to know how to plot a surface fitted to those points. I have found polynomial implementations of order 1 and 2 here: https://gist.github.com/amroamroamro/1db8d69b4b65e8bc66a6 but order 2 is too low.

My final goal is to obtain a volatility surface, such as this, extracted from this website:

In that blog post, python code is provided, but it is using NAG C propietary library. It says it is using Chebysev order 3 polynomial fit. I tried having a look at scipy's Chebyshev fit methods from numpy, but I haven't been able to understand how they work, and found no example.

Any help? I supposes this would be a piece of cake, but I am really lost here, since my math is a bit rusty.