I think you are on the right track in thinking about something like Zernike polynomials. As noted in the answer by jwimberly, these are an example of a system of orthogonal basis functions on a disk. I am not familiar with Zernike polynomials, but many other families of orthogonal functions (including Bessel functions) arise naturally in classical mathematical physics as eigenfunctions for certain partial differential equations (at the time of this writing, the animation at the top of that link even shows an example of a vibrating drum head).
Two questions come to my mind. First, if all you are after is the radial profile ($\theta$ averaged), then how much constraint on the spatial pattern do you need? Second, what types of variability occur in the spatio-temporal data?
In terms of the first question, there are two concerns that come to mind. Due to the polar coordinates, the support-area for each sensor has a trend with $r$. The second concern would be the possibility of aliasing, essentially a mis-alignment of your sensors relative to the phase of the pattern (to use a Fourier/Bessel analogy). Note that aliasing will likely be the primary uncertainty in constraining the peak temperatures (i.e. $T_{95}$).
In terms of this second question, data variability could actually help with any aliasing issues, essentially allowing any mis-alignment to average out over the different measurements. (Assuming no systematic bias ... but that would be a problem for any method, without e.g. a physical model to give more information).
So one possibility would be to define your spatial orthogonal functions purely at the sensor locations. These "Empirical Orthogonal Functions" could be computed via PCA on your spatiotemporal data matrix. (Possibly you could use some weighting to account for the variable sensor support areas, but given the uniform polar grid and target of radial averages, this may not be required.)
Note that if there is any physical modeling data available for "expected" variations in the temperature, available on a dense spatiotemporal computational grid, then the same PCA procedure could be applied to that data to derive orthogonal functions. (This would typically called "Proper Orthogonal Decomposition" in engineering, where it is used for model reduction, e.g. an expensive computational fluid dynamics model can be distilled for use in further design activities.)
A final comment, if you were to weight the sensor data by support area (i.e. polar cell size), this would be a type of diagonal covariance, in framework of GLS. (That would apply to your prediction problem more, although weighted PCA would be closely related.)
I hope this helps!
Update: Your new diagram of the sensor distribution changes things considerably in my view. If you want to estimate temperatures over the disk interior, you will need a much more informative prior than simply "set of orthogonal functions on the unit disk". There is just too little information in the sensor data.
If you indeed want to estimate the spatial temperature variation over the disk, the only reasonable way I can see would be to treat the problem as one of data assimilation. Here you would need to at least constrain the parametric form of the spatial distribution based on some physics-based considerations (these could be from simulations, or could be from related data in systems with similar dynamics).
I do not know your particular application, but if it is something like this, then I would imagine there is an extensive engineering literature that you could draw upon to choose appropriate prior constraints. (For that sort of detailed domain knowledge, this is probably not the best StackExchange site to ask on.)
