<p>You can create a list of RegionGeom objects for each point you want to plot and iteratively add them to your figure. Here is some code to do this</p>
<pre><code>import matplotlib.pyplot as plt
import astropy.units as u
from gammapy.maps import RegionGeom
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion, PointSkyRegion

ra = [248.3, 249.8, 250.2, 250.5, 250.5]
dec = [68.8, 67.7, 65.7, 72.2, 63.3]
points = list(zip(ra, dec))

center = SkyCoord(ra=270*u.deg, dec=66*u.deg)
angle = 10 * u.deg

fig = plt.figure()
region = RegionGeom(CircleSkyRegion(center, angle))
ax = region.plot_region()
for p in points:
    ra_p, dec_p = p
    region_p = RegionGeom(PointSkyRegion(center=SkyCoord(ra=ra_p*u.deg, dec=dec_p*u.deg)))
    region_p.plot_region(
    ax=ax,
    facecolor=&quot;black&quot;,
    edgecolor=&quot;black&quot;,
    kwargs_point={
        &quot;color&quot;: &quot;black&quot;,
        &quot;fillstyle&quot;: &quot;full&quot;,
        &quot;marker&quot;: &quot;o&quot;,
        &quot;markersize&quot;: 3,
    },
)
plt.grid()
plt.show()
</code></pre>
<p>Which produces this figure</p>
<p><img src="https://i.sstatic.net/yURdIQ0w.png" alt="image" /></p>

CC BY-SA 4.0
2024-09-18T22:25:35.137Z
*Welcome to Astronomy!* It's great a new user drops in and immediately answers a non-trivial question with a thorough answer.
uhoh
CC BY-SA 4.0
2024-09-19T09:08:23.487Z
thank you and welcome to Astronomy Stackexchange! Your solution works, but I hoped for something more elegant that would be beyond forcing a loop of points. I think such functionality exists. Nevertheless, I upvoted your answer :)
NeStack
CC BY-SA 4.0
2024-09-18T04:12:08.200Z
2024-09-18T04:18:08.340Z
2024-09-18T04:18:08.340Z
Austin Shen
<p>I am trying to do a scatter plot of a few data points on a region of the sky and then add a circle to the plot. Note, that I don't want a plot of the full sky, just of the region around the data points. The coordinates of the data points in degrees are:</p>
<pre><code>ra = np.array([248.3, 249.8, 250.2, 250.5, 250.5])
dec = np.array([68.8, 67.7, 65.7, 72.2, 63.3])
</code></pre>
<p>The centre of the circle should be at</p>
<pre><code>ra1 = 270
dec1= 66
</code></pre>
<p>It seems the packages for that (because we need to account for the curvature of the sky!) would be <code>astropy</code> and <code>regions</code>. But I just can't make it work. This post <a href="https://stackoverflow.com/a/77265601/5838180">here</a> is almost what I want to achieve, but it works with only one point and I don't see how to add more points. Thx for advice!</p>

CC BY-SA 4.0
2024-09-16T16:25:32.700Z
You probably just need matplotlib.  The fact that your data is RA/Dec coordinates doesn't really come in to play, just graph it as generic points.
Greg Miller
CC BY-SA 4.0
2024-09-16T17:00:31.277Z
My problem is that I want to check how many of my data points (that are more than the dummy points I listed) are within the circle and how many are outside. For that I need to draw a circle that is physically accurate projection on the sky and matplotlib doesn't seem to help me out for that task
NeStack
CC BY-SA 4.0
2024-09-17T02:24:25.760Z
Yes I think you are right; it certainly looks like you either need several of those packages, or you need to do the math yourself. "I  tried it and it didn't work" is usually not good enough to get a solution in Stack Exchange - it's necessary to show exactly what you tried and how it failed.
uhoh
CC BY-SA 4.0
2024-09-17T02:24:32.433Z
As for doing the math yourself, A "circle in the sky" is the intersection of a cone with the celestial sphere. It is not the same thing as a circle drawn on a cartesian plot of RA vs Dec. The circle would be formed by the collection of all unit vectors such that their dot product with a unit vector pointing at the center of the circle equals $\cos(\theta_0)$ where $\theta_0$ is the angular radius of your circle.
uhoh
CC BY-SA 4.0
2024-09-17T02:25:56.480Z
@GregMiller I disagree. See https://i.sstatic.net/kOWrz.png from the linked answer.
uhoh
CC BY-SA 4.0
2024-09-17T02:35:29.320Z
@uhoh I'm not going to disagree with you, but note that image is almost 30 degrees wide in RA and near the poles.  For the points provided, the distortion isn't really an issue, but for a more generalized solution it would be.
Greg Miller
<coordinate><python><astropy>
Code a scatter plot on a region of the sky with a circle
NeStack
# Code a scatter plot on a region of the sky with a circle

(Asked by: username_0 on 2024-09-16)

I am trying to do a scatter plot of a few data points on a region of the sky and then add a circle to the plot. Note, that I don't want a plot of the full sky, just of the region around the data points. The coordinates of the data points in degrees are:
ra = np.array([248.3, 249.8, 250.2, 250.5, 250.5])
dec = np.array([68.8, 67.7, 65.7, 72.2, 63.3])

The centre of the circle should be at
ra1 = 270
dec1= 66

It seems the packages for that (because we need to account for the curvature of the sky!) would be astropy and regions. But I just can't make it work. This post here is almost what I want to achieve, but it works with only one point and I don't see how to add more points. Thx for advice!

**Tags:** coordinate, python, astropy
> **Comments:**
> **username_1**: You probably just need matplotlib.  The fact that your data is RA/Dec coordinates doesn't really come in to play, just graph it as generic points.  
> **username_0**: My problem is that I want to check how many of my data points (that are more than the dummy points I listed) are within the circle and how many are outside. For that I need to draw a circle that is physically accurate projection on the sky and matplotlib doesn't seem to help me out for that task  
> **username_2**: Yes I think you are right; it certainly looks like you either need several of those packages, or you need to do the math yourself. "I  tried it and it didn't work" is usually not good enough to get a solution in Stack Exchange - it's necessary to show exactly what you tried and how it failed.  
> **username_2**: As for doing the math yourself, A "circle in the sky" is the intersection of a cone with the celestial sphere. It is not the same thing as a circle drawn on a cartesian plot of RA vs Dec. The circle would be formed by the collection of all unit vectors such that their dot product with a unit vector pointing at the center of the circle equals $\cos(\theta_0)$ where $\theta_0$ is the angular radius of your circle.  
> **username_2**: @GregMiller I disagree. See https://i.sstatic.net/kOWrz.png from the linked answer.  
> **username_1**: @uhoh I'm not going to disagree with you, but note that image is almost 30 degrees wide in RA and near the poles.  For the points provided, the distortion isn't really an issue, but for a more generalized solution it would be.  

---

## 1 answers

### Answer by username_3 on 2024-09-18. Score: 5

You can create a list of RegionGeom objects for each point you want to plot and iteratively add them to your figure. Here is some code to do this
import matplotlib.pyplot as plt
import astropy.units as u
from gammapy.maps import RegionGeom
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion, PointSkyRegion

ra = [248.3, 249.8, 250.2, 250.5, 250.5]
dec = [68.8, 67.7, 65.7, 72.2, 63.3]
points = list(zip(ra, dec))

center = SkyCoord(ra=270*u.deg, dec=66*u.deg)
angle = 10 * u.deg

fig = plt.figure()
region = RegionGeom(CircleSkyRegion(center, angle))
ax = region.plot_region()
for p in points:
    ra_p, dec_p = p
    region_p = RegionGeom(PointSkyRegion(center=SkyCoord(ra=ra_p*u.deg, dec=dec_p*u.deg)))
    region_p.plot_region(
    ax=ax,
    facecolor="black",
    edgecolor="black",
    kwargs_point={
        "color": "black",
        "fillstyle": "full",
        "marker": "o",
        "markersize": 3,
    },
)
plt.grid()
plt.show()

Which produces this figure

> **Comments:**
> **username_2**: *Welcome to Astronomy!* It's great a new user drops in and immediately answers a non-trivial question with a thorough answer.  
> **username_0**: thank you and welcome to Astronomy Stackexchange! Your solution works, but I hoped for something more elegant that would be beyond forcing a loop of points. I think such functionality exists. Nevertheless, I upvoted your answer :)  

---

