Are units of angle really dimensionless?

The answers are no and no. Being dimensionless or having the same dimension is a necessary condition for quantities to be "compatible", it is not a sufficient one. What one is trying to avoid is called category error. There is analogous situation in computer programming: one wishes to avoid putting values of some data type into places reserved for a different data type. But while having the same dimension is certainly required for values to belong to the same "data type", there is no reason why they can not be demarcated by many other categories in addition to that.

Newton meter is a unit of both torque and energy, and joules per kelvin of both entropy and heat capacity, but adding them is typically problematic. The same goes for adding proverbial apples and oranges measured in "dimensionless units" of counting numbers. Actually, the last example shows that the demarcation of categories depends on a context, if one only cares about apples and oranges as objects it might be ok to add them. Dimension is so prominent in physics because it is rarely meaningful to mix quantities of different dimensions, and there is a nice calculus (dimensional analysis) for keeping track of it. But it also makes sense to introduce additional categories to demarcate values of quantities like torque and energy, even if there may not be as nice a calculus for them.

As your own examples show it also makes sense to treat radians differently depending on context: take their category ("dimension") viz. steradians or counting numbers into account when deciding about addition, but disregard it when it comes to substitution into transcendental functions. Hertz is typically used to measure wave frequency, but because cycles and radians are officially dimensionless it shares dimension with the unit of angular velocity, radian per second, radians also make the only difference between amperes for electric current and ampere-turns for magnetomotive force. Similarly, dimensionless steradians are the only difference between lumen and candela, while luminous intensity and flux are often distinguished. So in those contexts it might also make sense to treat radians and steradians as "dimensional".

In fact, radians and steradians were in a class of their own as "supplementary units" of SI until 1995. That year the International Bureau on Weights and Measures (BIPM) decided that "ambiguous status of the supplementary units compromises the internal coherence of the SI", and reclassified them as "dimensionless derived units, the names and symbols of which may, but need not, be used in expressions for other SI derived units, as is convenient", thus eliminating the class of supplementary units. The desire to maintain a general rule that arguments of transcendental functions must be dimensionless might have played a role, but this shows that dimensional status is to a degree decided by convention rather than by fact. In the same vein, ampere was introduced as a new base unit into MKS system only in 1901, and incorporated into SI even later. As the name suggests, MKS originally made do with just meters, kilograms, and seconds as base units, this required fractional powers of meters and kilograms in the derived units of electric current however.

As @dmckee pointed out energy and torque can be distinguished as scalars and pseudo-scalars, meaning that under the orientation reversing transformations like reflections, the former keep their value while the latter switch sign. This brings up another categorization of quantities that plays a big role in physics, by transformation rules under coordinate changes. Among vectors there are "true" vectors (like velocity), covectors (like momentum), and pseudo-vectors (like angular momentum), in fact all tensor quantities are categorized by representations of orthogonal (in relativity Lorentz) group. This also comes with a nice calculus describing how tensor types combine under various operations (dot product, tensor product, wedge product, contractions, etc.). One reason for rewriting Maxwell's electrodynamics in terms of differential forms is to keep track of them. This becomes important when say the background metric is not Euclidean, because the identification of vectors and covectors depends on it. Different tensor types tend to have different dimensions anyway, but there are exceptions and the categorizations are clearly independent.

But even tensor type may not be enough. Before Joule's measurements of the mechanical equivalent of heat in 1840s the quantity of heat (measured in calories) and mechanical energy (measured in derived units) had two different dimensions. But even today one may wish to keep them in separate categories when studying a system where mechanical and thermal energy are approximately separately conserved, the same applies to Einstein's mass energy. This means that categorical boundaries are not set in stone, they may be erected or taken down both for practical expediency or due to a physical discovery.

Many historical peculiarities in the choice and development of units and unit systems are described in Klein's book The Science of Measurement.

Whenever I think about this problem I go back to one of Joel Spolsky's articles, "Making Wrong Code Look Wrong", which talks about Hungarian notation. Not only the useless kind of Hungarian notation, where variables are named in a way that describes their types (f_pos for a float, d_pos for a double, etc.) - this is "Systems Hungarian" in the article - but the original, practical kind, "Apps Hungarian", where the name of a variable reflects what kind of physical quantity it represents. x_pos and y_pos for horizontal and vertical position, for example. Or in an example that might be more relevant to your case, circ_length and rad_length for circumference and radius, respectively.

In Apps Hungarian, if you ever found yourself writing circ_length + rad_length, you would suspect that something was wrong because you shouldn't be adding circumference and radius. Even though they're dimensionally consistent, they're not compatible in some sense. You'd want to rewrite it as something like this:

circ_length + circ_from_rad(rad_length)

Unit systems are a (somewhat limited) physical equivalent of Apps Hungarian. We use different units to denote different variables that are not compatible and shouldn't be added together.

This might seem like a weird perspective on units at first. After all, it seems intuitively obvious that length and time are different in some way that, say, height and width aren't. But then relativity came along, and we learned that length and time actually are compatible, you just need to do the right conversion.

prime_t, prime_x = prime_from_normal(normal_t, normal_x)

or in other words,

$$\begin{pmatrix}ct' \\ x'\end{pmatrix} = \begin{pmatrix}\gamma & -\beta\gamma \\ -\beta\gamma & \gamma\end{pmatrix}\begin{pmatrix}ct \\ x\end{pmatrix}$$

Then quantum mechanics came along and we learned that energy and frequency are also kind of the same thing.

energy = energy_from_freq(frequency)

$$E = \hbar\omega$$

And then quantum gravity came along and somebody invented Planck units, and that goes all the way down the rabbit hole to the point where everything is just a number, and you can freely add masses and charges and forces.

Stay away from quantum gravity.

Anyway, if it's so easy to reduce the number of distinct units by showing how they can all be converted into each other, you can also reverse the process. You'd treat the unit systems we have as reduced versions of more complicated systems that distinguish between quantities we normally think of as being the same. Height and width, for example. You could have "height-meter" and "width-meter" as effectively separate units. Or, in your case, "circumference-meter" and "radius-meter", in which case you'd define $$1\ \mathrm{rad} = \frac{1\ \text{circumference-meter}}{1\ \text{radius-meter}}$$ In your system, this is not going to be the same as $\frac{1\ \text{height-meter}}{1\ \text{width-meter}}$. You can make it the same by turning all these units into meters, in which case you'd recover the metric system, but then you lose the extra contextual information provided by your unit system.

Here's a practical example: slope $m$ is defined as height over (horizontal) length, $\Delta y/\Delta x$, which means the units of slope are $$[m] = \frac{\text{height-meter}}{\text{length-meter}}$$ On the other hand, angle $\alpha$ is defined as circumference over radius, $$[\alpha] = \frac{\text{circumference-meter}}{\text{radius-meter}}$$ The relationship between these two is $$m = \tan\alpha$$ so in this augmented unit system, you know that the tangent takes in units of $\frac{\text{circumference-meter}}{\text{radius-meter}}$ and gives a result in units of $\frac{\text{height-meter}}{\text{length-meter}}$.

However, there is a problem. What if you're doing a calculation involving the transverse and recessional velocities of a star? (Perpendicular and parallel to your line of sight, respectively.) In that case, you still use the tangent function, but you might get a result in units of $\frac{\text{length-meter}}{\text{height-meter}}$. Strictly speaking, this probably means you should have a separate function that would produce this kind of output. In practice, you can take it too far. Giving everything separate units is often more trouble than it's worth.

So you'll have to strike a balance between the two extremes. Many people agree that having angles designated by a unit to preserve the contextual information that they are angles (not something else) is useful. You can meaningfully use that information with trig functions: a function like $\tan$ has to take as input an angle, or a "circular" ratio of lengths (circumference to radius or such), and give as output a "rectangular" ratio of lengths (height to width or vice-versa or some such thing). The radian might be a "fake" unit, sure, but in a way it's no more fake than a unit of velocity, or angular momentum, and it is useful information to maintain.

Here is an entertaining mathematical answer. (Or at least, I find it entertaining, anyway.)

Let us take seriously the idea that we can treat radians as a unit, and proceed from there. This means that when we write an expression like $\sin \theta$, the argument $\theta$ must have units of radians, whereas the result (I'll assume) is just a number without any units. Or to put it another, way, the expression $\sin^{-1} x$ has units of radians.

Now, since we can write $\sin\theta$ as $\frac{i}{2}(e^{-i\theta}-e^{i\theta})$, that means the argument to $e^\theta$ must also have units of radians. Or, to put it another way, any expression of the form $\ln x$ has units of radians, since the logarithm is the inverse of the exponential.

However, we can quickly get into trouble here, since the integral of the logarithm is given by $$ \int \ln x \,dx = x\ln x - x + C. $$ Assuming $x$ is a dimensionless number, the term $x\ln x$ has units of radians, but the $x$ term is dimensionless. So we've reached an inconsistency of the kind you were looking for, where a quantity in radians is added to a dimensionless quantity, and we never even had to include any physics at all.

This is not necessarily an unsolvable problem. benrg has a nice answer, in which he points out that you can write the solution to the above integral as $x\ln(x/e) + C = x\ln x - x\ln e + C$, with the point being that $\ln e$ is a constant with the value 1 but units of radians, so the units of the expression are radians overall. This seems to be consistent and I rather like it.

All this might be worth thinking about a bit more some time, on a rainy afternoon when there's nothing else to do. My point here was to show that it isn't straightforward to treat radians as a unit, even in the world of pure maths. benrg shows that it might nevertheless be possible to do it consistently, which I find interesting.