I have a spherical rendering, where the spherical coordinates $\phi$ and $\theta$ are represented by the x and y axis of the image (similar to how world maps work):
Now given a point on the image with the pixel coordinates p($\phi$, $\theta$) I want to calculate the average color of all the neighboring pixels on the sphere within a radius (radius could maybe more easily be defined as an angle, since we are thinking of a unit sphere). For this I need to sum up all of the pixels within the circle, however the circle on the sphere will not be a circle on my map. So the question is, how can I get all the pixel coordinates (= my $\phi$ and $\theta$ angles) within the circle?
Edit: Forgot to describe, how the rendering was done. With $\phi = [0, 2\pi]$ we divide $2\pi$ by the horizontal image resolution. $\theta = [0, \pi]$, so we divide $\pi$ by the vertical image resolution. Using a horizontal resolution that is twice the vertical resolution one pixel represents area on the sphere that is $(\frac{2 * \pi}{res_h}\times\frac{2 * \pi}{res_h})$.

I will stick to your notations $\phi$ for longitude, $\theta$ for latitude, although one finds rather often the inverse convention.
Given: a center $M_0$ (defined by spherical coordinates $(\phi_0,\theta_0)$) and a radius $R$ ($0<R<\pi/2$, measured as an arc on the unit sphere).
The set of points $M$ of the sphere which are interior to the spherical disk with center $M_0$ and radius $R$ is given by the following dot product constraint:
$$\vec{OM}.\vec{OM_0}>\cos(R)$$
which is equivalent, using classical spherical coordinates, to:
$$\begin{pmatrix}\cos(\phi)\cos(\theta)\\ \sin(\phi)\cos(\theta)\\ \sin(\theta)\end{pmatrix} .\begin{pmatrix}\cos(\phi_0)\cos(\theta_0)\\ \sin(\phi_0)\cos(\theta_0)\\ \sin(\theta_0)\end{pmatrix}>\cos(R)\tag{1}$$
This constraint can be written under the form:
$$\cos(\theta_0)\cos(\theta)\cos(\phi-\phi_0)+\sin(\theta_0)\sin(\theta)>\cos(R)\tag{2}$$
which is an implicit equation in $(\phi,\theta)$ depending upon three parameters $(\phi_0,\theta_0,R)$ that can be visualized with this Geogebra animation (play with the sliders !):
https://www.geogebra.org/calculator/eanc7njp
The top sliders $f$ and $t$ refer to the coordinates $\phi_0$ and $\theta_0$ resp. of center $M_0$.
Here are two examples:
Fig. 1: An "ordinary" circle centered in (0,0) with radius $\pi/4$ rendered as a kind of ellipse.
Fig. 2: A "limit case" image of a circle belonging to northern hemisphere, tangent to the equator, (almost) passing throughout North Pole, explaining the almost linear segment ranging approximately from $-\pi/2$ to $\pi/2$.
Remark 1: constraint (1) has been given in the same form by @blamocur.
Remark 2: formula (2) could have been obtained directly by using the spherical law of cosines..
Remark 3: Explicit equations of the circle can be found here.