First convert to polar coordinates $(x,y,t) \to (\rho,\phi,t) $. The equation becomes
$$ \frac{1}{c^2}u_t = u_{\rho\rho} + \frac{1}{\rho}u_{\rho} + \frac{1}{\rho^2}u_{\phi\phi} $$
Observe that $u = h$ is a steady-state solution satisfying the two boundary conditions, so we can write $u = h + v$, where $v$ also satisfies the heat equation but is homogeneous on the boundary
$$ \begin{align} v(r,\phi,t) &= 0 \\ v_\rho(R,\phi, t) &= 0 \\ v(\rho,\phi,0) &= H(\rho,\phi) - h \equiv f(\rho,\phi) \end{align} $$
Now use separation of variables. Let $v(\rho,\phi,t) = P(\rho)\Phi(\phi)T(t)$, then
$$ \frac{P''}{P} + \frac{P'}{\rho P} + \frac{\Phi''}{\rho^2 \Phi} = \frac{T'}{c^2 T} = -\lambda^2 $$
which gives $T(t) = e^{-c^2\lambda^2t}$. Next, separate again
$$ \frac{\rho^2P''}{P} + \frac{\rho P'}{P} + \lambda^2 \rho^2 = -\frac{\Phi''}{\Phi} = n^2 $$
Here, $n$ must be an integer to satisfy the periodicity condition $\Phi(\phi + 2\pi) = \Phi(\phi)$. Hence, $\Phi(\phi) = D\cos(n\phi) + E\sin(n\phi)$
For the radial function
$$ \rho^2 P'' + \rho P' + (\lambda^2\rho^2 - n^2)P = 0 $$
The substitution $z=\lambda\rho$ leads to the Bessel equation, thus we can write
$$ P(\rho) = A J_n(\lambda \rho) + B Y_n(\lambda \rho) $$
As stated in the comments, you have two separate problems due to the boundary condition at $\rho = r$. The solution on the smaller disk $\rho \in (0,r)$ is finite at the origin, so $B=0$ and
$$ P_{n,m}(\rho) = J_n(\lambda_{n,m} \rho) $$
where the eigenvalue $\lambda_{n,m}$ is such that $P(r) = J_n(\lambda_{n,m} r) = 0$. Here $m$ refer to the order of the zeroes of $J_n$, which can be found numerically.
The solution on the annulus $\rho \in (r,R)$ has different eigenvalues, so I'll call it $\mu$ to avoid confusion. The radial function needs to satisfy $P(r) = P'(R) = 0$ therefore
\begin{align}
A J_n(\mu r) + B Y_n(\mu r) &= 0 \\
A J_n'(\mu R) + B Y_n'(\mu R) &= 0
\end{align}
Eliminating either $A$ or $B$ leads to
$$ J_n(\mu r) Y_n'(\mu R) - Y_n(\mu r)J_n'(\mu R) = 0 $$
This is equation for $\mu$ that can be solved numerically. Once again, there are infinitely many solutions, so $\mu = \mu_{n,m}$. The derivatives of the Bessel function are given by the recurrence relations:
\begin{cases} J_0'(x) = -J_1(x) \\ J_n'(x) = \dfrac{J_{n-1}(x) - J_{n+1}(x)}{2} \end{cases}
Then the radial function is
$$ P_{n,m}(\rho) = Y_n(\mu_{n,m} r) J_n(\mu_{n,m} \rho) - J_n(\mu_{n,m} r) Y_n (\mu_{n,m} \rho) $$
We can assemble the general solution
$$ v(\rho,\phi,t) = \sum_{n=0}^\infty\sum_{m=0}^\infty P_{n,m}(\rho) \big[D_{n,m} \cos(n\phi) + E_{n,m}\sin(n\phi)\big]e^{-c^2\lambda_{n,m}^2 t} $$
Note that the function is piece-wise with different eigenvalues.
Matching the initial condition gives
$$ u(\rho,\phi,0) = f(\rho,\phi) = \sum_{n=0}^\infty\sum_{m=0}^\infty P_{n,m}(\rho) \big[D_{n,m} \cos(n\phi) + E_{n,m}\sin(n\phi)\big] $$
You can use orthogonality to find the constants. For example
$$ D_{n,m} = \dfrac{\displaystyle \int_0^{2\pi}\int_0^r f(\rho,\phi) J_n(\lambda_{n,m} \rho) \cos(n\phi) \ \rho\ d\rho\ d\phi}{\displaystyle \int_0^r \rho J_n^2(\lambda_{n,m}\rho)\ d\rho \int_0^{2\pi} \cos^2 (n\phi)\ d\phi} $$
