Indeed, the reshaping did not work. Always good to check shapes:
julia> r = reshape(x, :, 1, 1);
julia> size(r)
(13952, 1, 1)
The closest analog to numpy would probably be
const newaxis = [CartesianIndex()]
r = x[:, newaxis, newaxis, :]
θ = y[newaxis, :, newaxis, :]
ϕ = z[newaxis, newaxis, :, :]
X = vec(r .* sin.(θ) .* cos.(ϕ))
Alternatively, you can use something like Tullio.jl or TensorCast.jl and write it as
using TensorCast
@cast X[i⊗ j⊗ k⊗ b] := x[i, b] * sin(y[j, b]) * cos(z[k, b])
which closely matches your intended loop structure.