This can be realized by Multiwfn. After booting up Multiwfn and load a wavefunction file, enter the following commands
1  // Enter main function 1
d  // Decompose real space function at a point to orbital contributions
[Input X,Y,Z of a nucleus]
1  // Unit is Bohr (also you can choose 2 to select Angstrom)
5  // Spin density
[Press ENTER button]  // Contribution from 10 orbitals having largest contributions will be shown

Then you can see

Contribution from orbital    20 (occ= 1.000000): -0.13833944E+03 a.u.
Contribution from orbital     1 (occ= 1.000000):  0.13820669E+03 a.u.
Contribution from orbital     2 (occ= 1.000000):  0.63424863E+01 a.u.
Contribution from orbital    21 (occ= 1.000000): -0.60051753E+01 a.u.
Contribution from orbital     6 (occ= 1.000000):  0.12114241E+01 a.u.
Contribution from orbital    23 (occ= 1.000000): -0.57635741E+00 a.u.
Contribution from orbital     5 (occ= 1.000000):  0.57405887E+00 a.u.
Contribution from orbital    22 (occ= 1.000000): -0.67708475E-35 a.u.
Contribution from orbital     4 (occ= 1.000000):  0.95628644E-65 a.u.
Contribution from orbital     3 (occ= 1.000000):  0.00000000E+00 a.u.
Sum of above values:  0.14136869E+01 a.u.
Exact value:  0.14136869E+01 a.u.

The "exact value" is actual spin density at the given point, "Sum of above values" is sum of the printed terms.

Please also note that the index of beta orbitals is after that of alpha orbitals. After loading the file, you can see such as following information from screen

Orbitals from      1 to    19 are alpha, from     1 to     6 are occupied
Orbitals from     20 to    38 are beta,  from    20 to    23 are occupied

Therefore e.g. orbitals 20 and 21 are the first and second beta orbitals, respectively.