Indeed.

One way to fix the problem is to sample uniformly not on the latitude x longitude rectangle but the sin (latitude) x longitude rectangle.

The reason this works is because the area of a infinitesimal lat long patch on the sphere is dlong x lat x cosine (lat). Now, if we sample on the long x sin(lat) rectangle, an infinitesimal rectangle also has area dlong x dlat x d/dlat sin(lat) = dlong x dlat cos (lat).

Unfortunately, these simple fixes do not generalize to arbitrary dimensions. For that those that exploit rotational symmetry of L2 norm works best.