Random sample of points on a region of the surface of a unit sphere

I need to generate a random sample of points distributed on a region of the surface of a sphere. This answer gives an elegant recipe to sample the entire surface of the sphere:

def sample_spherical(npoints, ndim=3):
    vec = np.random.randn(ndim, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

(vec returns the (x, y, z) coordinates) which results in:

I need to restrict the sampling to a rectangular region. This region is defined by center values and lengths, for example:

• center coordinates, eg: (ra_0, dec_0) = (34, 67)
• side lengths, eg: (ra_l, dec_l) = (2, 1)

where (ra_0, dec_0) are coordinates in the equatorial system, and the lengths are in degrees.

I could use the above function in the following way: call it using a large npoints, transform the returned (x, y, z) values into the equatorial system, reject those outside the defined region, and repeat until the desired number of samples is reached. I feel that there must be a more elegant way of achieving this though.