Efficient least squares fit of 2D polynomials to data?

I have some data on a regular 2D grid that I would like to decompose into an orthogonal basis function set. What is the most efficient way to do the fit?

Here's some sample code:

import numpy as np

sze = (128,128)
x, y = np.linspace(-1, 1, sze[1]), np.linspace(-1, 1, sze[0])
xv, yv = np.meshgrid(x, y)
rho = np.sqrt(xv**2 + yv**2)
phi = np.arctan2(yv, xv)
data = (rho * np.cos(phi))
data[rho > 1] = 0
# now have some synthetic data -- need to fit to a basis function set that has single-index members:
# 0: 1
# 1: rho * cos(phi)
# 2: rho * sin(phi)
# ...
# 24: 70 rho^8 - 140 rho^6 + 90 rho^4 - 20 rho^2 + 1