Zhen Lin explained the theory here. I try to recite this in the theory section and then show ways to plot things with different tools in the Getting Our Hands Dirty! section.
Theory
The square form i.e. the Jacobian is required for the change in coordinates because you need to know change in every direction. After this calculation, you need to calculate the determinant and take the square root as below
$J = \begin{pmatrix} r \cos \theta \cos \phi & -r \sin \theta \sin \phi \\ r \sin \theta \cos \phi & r \cos \theta \sin \phi \\ -r \sin \phi & 0 \end{pmatrix}$
for which the general formula
$A = \int_S \sqrt{\det g} \, du \, dv.$
I really recommend Zhen Lin's marvelous answer!
P.s. The only case when your first equation is correct is when $r=1$ and $\sin(\phi)=1$ for the integration -- probably not the integral you are trying to do! I use the same notation as Zhen Lin, Latitude is $\psi \in [0,\pi]$ and azimuth is $\theta \in[0,2\pi]$.
ALERT!
Please, notice that people use two different conventions for integration: Zhen Lin uses reading outer-most integral things first and you use reading things in parallel (I hope you know the difference or you are very confused!).
Now let's get our hands dirty!
Now how should we alter the parametrization to make it look like the integral? I don't know, thinking. Anyway, the below examples were provided by co-operation of people around Stackexchange such as Anon Mouse and eindandi. If you find them useful, give them some upvotes -- not me.
Mathematica

Manipulate[ ParametricPlot3D[{Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]}, {theta, 0, t}, {phi, 0, p}, PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}], {t, 0.1, Pi}, {p, 0.1, 2 Pi}]
Matlab

[f,t] = meshgrid(linspace(0,2*pi,361),linspace(0,pi,361)); x = sin(t)*cos(f); y = sin(t)*sin(f); z = cos(t); surf(x,y,z)