for each point $p$ in theory we have a 3D function that tells us the orientation in a given direction
You have missunderstood this. $D$ is a Normal Distribution Function (or short NDF), so it doesn't really give you a single normal, but a distribution. In a (specular) BRDF you are always using the normal that is the half vector between the incoming and outgoing light, since via your theory, every microfacet is a perfect mirror and thus every microfacet reflects light exactly along the (micro) surface and that one only.
You also must not confuse roughness with microfacet normal. The roughness is a more or less arbitrary value (and a scalar, so a single one at that). The roughness is being used differently, depending on the overall BRDF (e.g. roughness in GGX(/Trowbridge-Reitz) is different from the Roughness in Oren-Nayar) and therefore can have different ranges. Still, if you look at one specific NDF, you will see, how the roughness is used.
$D(\omega_m) = \frac{\alpha^2}{\pi((\omega_n \cdot \omega_m)^2 (\alpha^2-1)+1)^2}$
with $\alpha = roughness^2$ (a common remapping), $\omega_m$ is the microfacet normal, $\omega_n$ is the geometry normal.
You can see here, that the roughness is being treated like a single value. The higher the roughness is, the more random your microfacets (or normals thereof) are distributed. The more that happens, the less your surface will reflect light concentrated into the same direction (i.e. the less highlights you will have). Thus, your material is more diffuse.
Since now your roughness is a single value, the map makes sense to have only values in $\left[0, 1\right]$ and therefore is a gray scale image.