If I correctly understand, here is how I would do in R.

library(uniformly)
library(pracma)
library(rgl)

# sample points on a spherical cap
points_on_cap1 <- runif_on_sphericalCap(300, r = 2, h = 0.5)
# convert to spherical coordinates
sphcoords1 <- cart2sph(points_on_cap1)

# sample points on a spherical cap
points_on_cap2 <- runif_on_sphericalCap(300, r = 2, h = 0.5)
# rotate them, because this is the same spherical cap as before
points_on_cap2 <- rotate3d(points_on_cap2, 3*pi/4, 1, 1, 1)
# convert to spherical coordinates
sphcoords2 <- cart2sph(points_on_cap2)

# 3D plot
spheres3d(0, 0, 0, radius = 2, alpha = 0.5, color = "yellow")
points3d(points_on_cap1, color = "blue")
points3d(points_on_cap2, color = "red")

# 2D plot (of the spherical coordinates)
plot(
  sphcoords1[, 1:2], xlim = c(-pi, pi), ylim = c(-pi/2, pi/2), 
  pch = 19, col = "blue"
)
points(sphcoords2[, 1:2], pch = 19, col = "red")

Do I understand?

Here is the function runif_on_sphericalCap:

function(n, r = 1, h){
  stopifnot(h > 0, h < 2*r)
  xy <- runif_in_sphere(n, 2L, 1)
  k <- h * apply(xy, 1L, crossprod)
  s <- sqrt(h * (2*r - k))
  cbind(s*xy, r-k)
}

It always samples on a spherical cap with symmetry axis joining the center of the sphere to the North pole. That is why I do a rotation, to get another spherical cap.

Say me if I understand and I'll try to help you to convert the code to Julia.

EDIT: Julia code

using Random, Distributions, LinearAlgebra

function runif_in_sphere(n::I, d::I, r::R) where {I<:Integer, R<:Number}
  G = Normal()
  sims = rand(G, n, d)
  norms = map(norm, eachrow(sims))
  u = rand(n) .^ (1/d)
  return r .* u .* broadcast(*, 1 ./ norms, sims)
end

function runif_on_sphericalCap(n::I, r::Number, h::Number) where {I<:Integer}
  if h <= 0 || h >= 2*r
    error("")
  end
  xy = runif_in_sphere(n, 2, 1.0)
  k = h .* map(x -> dot(x,x), eachrow(xy))
  s = sqrt.(h .* (2*r .- k))
  return hcat(broadcast(*, s, xy), r .- k)
end