As for asking about whether the symmetries of this equation would help you solve it, here are a few remarks that you may (or may not) find useful:
I assume that you want to consider what are usually called the 'point symmetries', i.e., the transformations of the $r\theta$-plane that carry (graphs of) solutions of the equation $$ r'(\theta)^2 + r(\theta)^2 = \theta^2 \tag{1} $$ into (graphs of) solutions of this equation. Now, by definition, the graphs of solutions of (1) are the null curves of the quadratic differential $$ g = dr^2 - (\theta^2{-}r^2)\ d\theta^2 $$ and, of course, this quadratic differential is positive definite in the open domain consisting of the two wedges where $r^2 > \theta^2$, so there are no null curves in this region. Meanwhile, in the open region $U$ described by the inequality $r^2 <\theta^2$, the quadratic form $g$ factors as $$ g = \bigl(dr - \sqrt{\theta^2{-}r^2}\ d\theta \bigr)\bigl(dr + \sqrt{\theta^2{-}r^2}\ d\theta \bigr), $$ So this region is foliated by two transverse foliations of $g$-null curves, namely, the foliation $F_+$ consisting of the solutions to $$ \frac{dr}{d\theta} = \sqrt{\theta^2{-}r^2} $$ and the foliation $F_-$ consisting of the solutions to $$ \frac{dr}{d\theta} = -\sqrt{\theta^2{-}r^2} $$ Obviously, the transformation $(\theta,r)\mapsto(\theta,-r)$ exchanges these two foliations, while the transformation $(\theta,r)\mapsto(-\theta,-r)$ preserves them (and exchanges the two components of $U$). Thus, it's enough to describe $F_+$ in the region of $U$ where $\theta>0$.
Each of the leaves of $F_+$ in the region of $U$ where $\theta>0$ consists of a graph of an increasing solution of (1) that is defined on an interval $(|a|,\infty)$ and that satisfies $r\bigl(|a|\bigr) = a$ and $r'\bigl(|a|\bigr) = 0$ while $r'(t)>0$ for all $t>|a|$. When $a\not=0$, one has a parametric expression for $r$ near $\theta=|a|$ of the form $$ r = a + \frac{|a|^{1/2}}{3}\ s^3 -\frac{a}{12}\ s^4 + \frac{2a^2+9}{360|a|^{1/2}}\ s^5 + \cdots\qquad\text{and}\qquad \theta = |a| + \frac12\ s^2 $$ when $s>0$. (When $s < 0$ this is a leaf of $F_-$.) It is easy to see that each leaf of $F_+$ in the quadrant where $\theta >0$ is not bounded above by any fixed $r$-value and hence each leaf of $F_-$ in this region is not bounded below by any fixed $r$-value.
Thus, if a point $(\theta,r)$ satisfies $\theta>0$ and $\theta^2-r^2>0$, then there will be unique numbers $x(\theta,r) < y(\theta,r)$ such that the leaf of $F_+$ that passes through $(\theta,r)$ also passes through $\bigl(|x(\theta,r)|,x(\theta,r)\bigr)$ while the leaf of $F_-$ that passes through $(\theta,r)$ also passes through $\bigl(|y(\theta,r)|,y(\theta,r)\bigr)$.
The mapping $(\theta,r)\mapsto \bigl(x(\theta,r),y(\theta,r)\bigr)$ carries this quadrant in the $\theta r$-plane onto the half-plane $y>x$ in the $xy$-plane and carries the $F_+$ foliation to the $x=const$ foliation and the $F_-$ foliation into the $y=const$ foliation.
Now the point symmetry group is clear: It is the set of transformations in this open quadrant of the $\theta r$-plane that gets carried by this mapping into the transformations of the half-plane $y>x$ that are of the form $(x,y)\mapsto \bigl(f(x),f(y)\bigr)$ where $f$ is a strictly increasing diffeomorphism of the line.
However, I don't think that you will be able to find this symmetry group explicitly without integrating the equations, which, of course, you don't want to have to do because getting help in finding the solutions is why you are interested in the symmetries in the first place.
In any case, it is clear that there are two (and only two) solutions that satisfy $r(0)=0$. Each is the negative of the other, and one can easily find the Taylor series for the solution whose graph lies in $F_+$ when $\theta>0$ and $F_-$ when $\theta<0$. It is an even function of $\theta$ and satisfies $$ r(\theta) = \frac12\theta^2 - \frac{1}{32}\theta^4 + \frac{1}{768}\theta^6 - \frac{1}{49152}\theta^8 + \cdots $$ as described above in my comment.
Now, because the equation (1) is singular, it is not obvious that this formal power series converges; in fact, one needs a little theorem to prove this. The best way to get to this is to assume a solution in the form $$ r(\theta) = \theta^2\bigl(\tfrac12 + u(\theta^2)\bigr) $$ where $u(0)=0$. Then one finds that, setting $\theta^2=x$, one has an equation for $u(x)$ in the form $$ xu'(x) = \tfrac12\bigl(1 - x\bigl(\tfrac12 + u(x)\bigr)^2\bigr)^{1/2}-\tfrac12 - u(x) = H\bigl(x,u(x)\bigr) $$ where $H$ is an analytic function near $(0,0)$ that satisfies $H(0,0)=0$ and $H_u(0,0)=-1$ (in particular, $H_u(0,0)$ is not a positive integer). As such, this is in the standard Briot-Bouquet form, and thus one knows that there exists a unique analytic solution $u(x)$ in a neighborhood of $x=0$ that satisfies $u(0)=0$. (See R. Gérard, Étude locale des équations différentielles de la forme $xy'=f(x,y)$ au voisinage de $x=0$, J. Fac. Sci. Univ. Tokyo 36 (1989) and the references cited therein.) This proves the convergence of the formal series for $r(\theta)$.
With a little more work, one can show that this solution $r(\theta)$ exists for $-\infty<\theta<\infty$.
