@LightPhoton The explanation given in the book sweeps some important methodology under the rug, this methodology being a standard way to calculate the number of available quantum states of a system - not only for an ideal gas but also for black-body radiation (this gives you the famous Planck formula for its energy density) and virtually any system composed of weakly interacting parts. This standard calculation involves, as mentioned in the book, closing a system in a rectangular box, however at a later stage of the calculation you indeed introduce auxiliary variables which you deal with by "closing" them inside... a spherical box

. I suppose that your confusion caused by the book using seemingly "rectangular" as well as "spherical" boxes in the "momentum" as well as "coordinate" spaces is precisely due to it not showing the technique I'm speaking about explicitly, instead arriving at the (correct!) results as quickly as possible to move on to different topics and applications.
The explicit methodology is the following. Consider a large rectangular box with sides of some length ##L## and volume ##L^3##. The particles inside the box (your ideal gas) can occupy any positions ##\mathbf{r}=(x,y,z)##, while the momentum eigenstates of these particles are described by the usual exponential wave functions
$$
e^{\frac{i}{\hbar}\mathbf{p}\cdot\mathbf{r}} = e^{\frac{i}{\hbar}p_x x}e^{\frac{i}{\hbar}p_y y}e^{\frac{i}{\hbar}p_z z} \rm{.}
$$
Now, since you enclosed your system in a rectangular box, the system is no longer translationally-invariant (because of the presence of the walls). The principles of mechanics then tell you that actually the system's momentum is not conserved, so it is rather awkward to work with momentum vectors ##\mathbf{p}=(p_x,p_y,p_z)## as they will not actually be well-defined quantities. To remedy this - and this is one of the important steps - you enforce the so-called
periodic boundary conditions, namely you assume that when the particle's position coincides with one of the walls of the box, then this particle "leaves" the box and another identical particle "enters" the box through a wall directly opposite to that through which the "original" particle left. Basically, you artificially impose the translation-invariance back on your system, so that you may still speak about momenta ##\mathbf{p}## while dealing with a finite (but large) box. Since the lengths of all sides of the box are equal, the periodic boundary conditions tell you that the state of the particle located at ##(x,y,z)## is the same as if it was located at ##(x+L, y+L, z+L)##. Thus
$$
e^{\frac{i}{\hbar}p_x (x+L)}e^{\frac{i}{\hbar}p_y (y+L)}e^{\frac{i}{\hbar}p_z (z+L)} = e^{\frac{i}{\hbar}p_x x}e^{\frac{i}{\hbar}p_y y}e^{\frac{i}{\hbar}p_z z}
$$
which shows that you must have
$$
p_x = \frac{2\pi\hbar}{L}n_x, \quad p_y = \frac{2\pi\hbar}{L}n_y, \quad p_z = \frac{2\pi\hbar}{L}n_z
$$
where ##n_x##, ##n_y## and ##n_z## are integers. Now, since the momentum
vectors ##\mathbf{p}## can have arbitrary direction, the three integers ##n_i## can be either positive or negative numbers, and based on the above the momentum has the form
$$
\mathbf{p} = (p_x, p_y, p_z) = \frac{2\pi\hbar}{L}(n_x,n_y,n_z) = \frac{2\pi\hbar}{L}(\pm|n_x|,\pm|n_y|,\pm|n_z|) \rm{.}
$$
So you obtain that to every triple ##(|n_x|,|n_y|,|n_z|)## of
positive integers there corrsponds ##2^3=8## momentum vectors ##\mathbf{p}##.
Now comes the part I mentioned eariler, with the auxiliary variables which will be "enclosed" in a spherical box. Namely, you see that the problem of counting the number of availalbe states for given maximum value of momentum ##\mathbf{p}_\text{max}## can be reduced, in this case, to counting the number of the integer triples ##(|n_x|,|n_y|,|n_z|)## for given maximum value of these triples. You do this counting by drawing a three-dimensional grid of points in the coordinate system whose three axes are labeled with ##|n_x|##, ##|n_y|## and ##|n_z|## (instead of, say, ##x##, ##y## and ##z##). Now, if there are sufficiently many such points then you can draw a sphere centered in the origin of your system of axes and this sphere will contain many points ##(|n_x|,|n_y|,|n_z|)##. If you choose some maximum values of ##|n_x|##, ##|n_y|## and ##|n_z|## that you wish to consider, then the number of distinct triples of integers that you need to count will be given by the number of points enclosed by the sphere of radius ##\sqrt{n_x^2 + n_y^2 + n_z^2}##. If the sphere is large enough, then this number is given by one eighth (1/8) of the volume of your sphere. The 1/8 factor comes from the fact that you consider only positive integers ##|n_x|##, ##|n_y|## and ##|n_z|##, so that only one octant of the sphere corresponding to positive semi-axes of your coordinate system contribute to the counting.
Thus the number of distinct triples of positive integers ##(|n_x|, |n_y|, |n_z|)## enclosed by the sphere is
$$
\frac{1}{8} \cdot \frac{4}{3}\pi\left(\sqrt{n_x^2 + n_y^2 + n_z^2}\right)^3
$$
Now, remember that earlier we showed that to each triple of positive integers ##(|n_x|, |n_y|, |n_z|)## there corresponds ##2^3=8## momentum vectors ##\mathbf{p}##. Since the above expression is the number of such integer triples, to obtain the corresponding number of momentum vectors ##\mathbf{p}## you need to multiply the above equation by 8. You also know that
$$
p^2 = \mathbf{p}\cdot\mathbf{p} = \left(\frac{2\pi\hbar}{L}\right)^2\left(n_x^2 + n_y^2 + n_z^2\right)
$$
and using this you can finally obtain that the desired number of momentum vectors whose magnitude is at most ##p## is given by
$$
N_p = \frac{4\pi}{3} \frac{L^3}{(2\pi\hbar)^3} p^3 \rm{.}
$$
Now, this is the total number of momenta corresponding to the number of integer triples inside the sphere in the "##n_i##-space". To obtain the number of momentum vectors ##\text{d}N_p## whose magnitude lies between ##p## and ##p + \text{d}p## (as is usually done) you differentiate the above expression with respect to ##p##, obtaining:
$$
\text{d}N_p = \frac{\partial N_p}{\partial p}\text{d}p = \frac{L^3}{(2\pi\hbar)^3} 4\pi p^2 \text{d}p \rm{.}
$$
You may recognize the ##4\pi## factor above as just the total solid angle, and in fact
$$
4\pi p^2 \text{d}p = \int_0^\pi \int_0^{2\pi} p^2 \sin\Theta_p \text{d}p \,\text{d}\Theta_p \,\text{d}\phi_p \rm{,}
$$
where ##p##, ##\Theta_p## and ##\phi_p## are the spherical components of the vector ##\mathbf{p}##. Thus the number ##\text{d}N_p## obtained above is just the expression
$$
\text{d}\eta_\mathbf{p} = \frac{L^3}{(2\pi\hbar)^3} p^2 \sin\Theta_p \text{d}p \,\text{d}\Theta_p \,\text{d}\phi_p = \frac{L^3}{(2\pi\hbar)^3} \text{d}^3p
$$
integrated over all angles.
This is the final answer. The number of the available momentum states for which the components of ##\mathbf{p}## lie in the intervals ##(p_x, p_x+\text{d}p_x)##, ##(p_y, p_y+\text{d}p_y)## and ##(p_z, p_z+\text{d}p_z)## is ##\frac{L^3}{(2\pi\hbar)^3} \text{d}^3p##. From this you get that the
density of states, that is the number of states divided by the volume of the box ##L^3##, is ##\frac{\text{d}^3p}{(2\pi\hbar)^3}##.
Of course since the particles in the box can take any position ##\mathbf{r}##, it is obvious that the number of states for which the coordinates of the particles lie between ##(x, x+\text{d}x)##, ##(y, y+\text{d}y)## and ##(z, z+\text{d}z)## is just ##\text{d}^3r = \text{d}x \,\text{d}y \,\text{d}z##. Combining this with the result we derived above, you simply get that the density of states of your gas (the number of available states per unit volume of the box - see? at the end the presence of the box doesn't matter

) is
$$
\frac{\text{d}^3r \,\,\text{d}^3p}{(2\pi\hbar)^3} = \frac{\text{d}^3r \,\,\text{d}^3p}{h^3}
$$
where I just used the fact that ##2\pi\hbar = h##, to be consistent with the notation used by the book.
As a bonus I will add a side comment explaining why the above result is remarkable (at least to me). Observe namely that the formula derived above is dimensionless, because of the presence of the Planck constant ##h## in the denominator. As
@anuttarasammyak already pointed out, the "bare" product ##\text{d}x \,\text{d}p_x## of the position-momentum variables has a dimension, namely those of action ##J\cdot s## (Joule times second). Now, in statistical mechanics one of the most important concepts is that of the entropy, which is given as the logarithm of the number of states available to the system. In classical mechanics, where there is no Planck constant, the usual phase-space formulation gives you the density of states precisely as this "bare" product ##\text{d}x \,\text{d}p_x##. Therefore the entropy ##S = \ln{\left(\Delta x \,\Delta p_x\right)}## is defined only up to an additive constant, which depends on the choice of units one adopts for energy and time (or position and momentum): ##S = \ln{\left(\alpha\left[\Delta x \,\Delta p_x\right]\right)} = \ln{\left(\Delta x \,\Delta p_x\right)} + \ln{\alpha}##. Because of this, in classical mechanics only entropy
differences can have any physical meaning, as these additional constant ##\ln{\alpha}## terms cancel each other upon taking the difference. Observe, however, that quantum mechanics makes possible - by introducing the Planck constant - to define the entropy in an absolute sense, ##S = \ln{\left(\frac{\Delta x \,\Delta p_x}{h}\right)}## as in this case the argument of the logarithm does not depend on the choice of units. Just a little neat observation.