We present an algorithm for automatic integration over an N-dimensional sphere. The quadrature formula is obtained by using a trapezoidal rule after a non-linear transformation, and allows to deal with integrand singularities on the surface or in the centre of the sphere. At the basis of the theoretical development lie the construction and the selection of suitable transformations. The algorithm is cast into an automatic integration program coded as a Standard Fortran sub-routine.