It isn't uniform in colatitude and longitude. (Sorry, I have to use those terms because I use spherical coordinates with theta and phi the opposite way round to you. Colatitude is angle from the north pole. Longitude is angle from one meridion.)
A small increment of surface on the sphere r=1 is
sin(colatitude) d(colatitude) d(longitude)
or
d(cos(colatitude)) d(longitude)
or
d(z) d(longitude)
Ah, I see what @Thomas Huxhorn means now  it is uniform in z and longitude, with z ranging from 1 to 1 and longitude from 0 to 2.pi, so giving a total surface area of 4 pi (which is correct).
So, generate z uniform in [1,1]
and longitude uniform in [0, 2 pi]
You can recover colatitude from cos(colatitude)=z, which is where the acos() function is coming from.
Grrr. I code in Fortran most of the time but it took me too long to decipher that. Note that Fortran's standard randomnumber generation is uniform on [0,1) (although it's useful that you can get an entire arrayful at once, rather than repeated calls).
Here's an attempt in C++.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

#include <iostream>
#include <iomanip>
#include <random>
#include <chrono>
#include <cmath>
using namespace std;
const double PI = 3.141592653589793238;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
mt19937 gen( seed );
uniform_real_distribution<double> dist( 0.0, 1.0 );
//====================================================================
// Return random point on the unit sphere r=1
// Note that it returns both the cartesian coordinates x, y, z
// and the spherical angles colatitude and longitude (in radians)
//====================================================================
void getPoint( double &x, double &y, double &z, double &colatitude, double &longitude )
{
z = 1.0 + 2.0 * dist( gen ); // uniformly distributed z
longitude = 2.0 * PI * dist( gen ); // uniformly distributed longitude
colatitude = acos( z );
double rh = sqrt( 1.0  z * z ); // horizontal radius (= r sin(colatitude), with r=1 )
x = rh * cos( longitude );
y = rh * sin( longitude );
}
int main()
{
#define SP << setw(10) << setprecision(4) << fixed <<
const int N = 10;
double x, y, z, colatitude, longitude;
cout SP "x" SP "y" SP "z" << '\n';
for ( int i = 0; i < N; i++ )
{
getPoint( x, y, z, colatitude, longitude );
cout SP x SP y SP z << '\n';
}
}
 
x y z
0.8246 0.1450 0.5469
0.2953 0.0332 0.9548
0.9162 0.0761 0.3934
0.2246 0.5299 0.8178
0.6054 0.2804 0.7449
0.5632 0.6979 0.4424
0.3306 0.7446 0.5799
0.3448 0.7311 0.5888
0.6028 0.0653 0.7952
0.0447 0.6338 0.7722 