As lightxbulb mentioned in the comments, one confusing thing here is that two different functions are being called $p$. Kind of like function overloading in C++, the $p(\omega)$ and $p(\theta, \phi)$ are two different functions, in the sense of having different domains and different dependence on their inputs, though in some sense they abstractly represent the same distribution.
The other thing that's confusing here is that the values returned by these functions are not just numbers. Since these are probability density functions, their values carry units of inverse area. Moreover, the two functions are expressing density with respect to different areas: $p(\omega)$ is probability per area on the unit sphere, while $p(\theta, \phi)$ is probability per unit area on the $\theta$-$\phi$ plane—an abstract 2D plane formed by $\theta$ and $\phi$ as Cartesian coordinates (just like the $x$-$y$ plane). These area measures are related to each other, but not the same: as you observed, $d\omega = \sin \theta \, d\theta \, d\phi$.
Notation like $d\omega$ and $d\theta \, d\phi$ refers to these area measures. As mentioned in the comments, this can be mathematically formalized as differential forms (here's an introduction I found useful).
Because $p(\omega)$ is a probability density per area on the unit sphere, if $d\omega$ is a "small" area on the unit sphere around the point $\omega$, then $p(\omega) \, d\omega$ works out to be a probability—not a density, but just a plain old probability number. In the limit of $d\omega$ being infinitesimal, $p(\omega) \, d\omega$ is the amount of probability enclosed within the area $d\omega$.
Similarly $p(\theta, \phi) \, d\theta \, d\phi$ is the amount of probability enclosed within a "small" area of the $\theta$-$\phi$ plane. The statement that
p(\omega) \, d\omega = p(\theta, \phi) \, d\theta \, d\phi
is intended to mean the following: the amount of probability in a small area of the unit sphere is equal to the amount of probability in the corresponding small area of the $\theta$-$\phi$ plane. In other words, probability is conserved by mapping between these two domains. This is the condition that makes these two functions represent "the same" distribution, abstractly, despite being parameterized with different variables.
The condition needs to be written that way because the probability densities aren't unitless numbers. If we had some function whose values were just plain numbers without units attached, and we wanted to reparameterize it with different variables, the condition would just be that $f(\omega) = f(\theta, \phi)$: the function should keep the same value at "the same" point when we change variables. But with the probability densities, we don't want the same value of the function, we want the same amount of probability in each part of the domain.
Using this "conservation of probability" condition $p(\omega) \, d\omega = p(\theta, \phi) \, d\theta \, d\phi$ together with the relationship between the two area measures, $d\omega = \sin \theta \, d\theta \, d\phi$, we can derive:
p(\theta, \phi) = p(\omega) \, \sin \theta
which shows that to express the a spherical distribution in terms of $\theta, \phi$, you have to not only change variables but also multiply the function by $\sin \theta$ to conserve the probability. That $\sin \theta$ factor is known as the "inverse Jacobian", and it accounts for how the transformation stretches the area in different parts of the domain.