<p>As long as we can live with a few tens of kilometers accuracy, we can treat the Earth's surface as a sphere.</p>
<p>The subpoint of the object (the point on Earth at which the object appears at the zenith, overhead, or the central point on Earth seen from the view of the object if they could see Earth clearly) is <a href="https://en.wikipedia.org/wiki/Hour_angle#Relation_with_the_right_ascension" rel="nofollow noreferrer">simply given by</a></p>
<p><span class="math-container">$$\text{lat} = \text{Dec}$$</span></p>
<p><span class="math-container">$$\text{lon} = \text{R.A.} - \text{GST}$$</span></p>
<p>where GST is <a href="https://en.wikipedia.org/wiki/Sidereal_time" rel="nofollow noreferrer">Greenwich Sidereal time</a></p>
<p>To make your plot we draw a great circle around the subpoint, then make a simple <a href="https://en.wikipedia.org/wiki/Equirectangular_projection" rel="nofollow noreferrer">Equirectangular projection</a> or lat/lon projection.</p>
<p>In Python for example, calculate the latitude and longitude of the subpoint,</p>
<pre><code>latsub = np.radians(dec)
lonsub = np.radians(RA - GST) # where these have already been converted to degrees
</code></pre>
<p>Pick a bunch of points to make a circle</p>
<pre><code>theta = np.linspace(0, 2*np.pi, 10001)
</code></pre>
<p>Draw a unit circle around Earth's equator, assuming the subpoint is at the north pole (dec = +90 degrees).</p>
<pre><code>x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
</code></pre>
<p>Tilt it by 90 - dec</p>
<pre><code>tilt = np.pi/2 - np.radians(dec)
x, z = [x * f(tilt) for f in (np.cos, np.sin)]
</code></pre>
<p>calculate the latitude and longitude of the tilted circle by shifting the longitude (subtracting the longitude of the subpoint)</p>
<pre><code>loncirc = np.arctan2(y, x) - lonsub
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))
</code></pre>
<p>Put it all together and convert units properly, and you get something like this. I did it for the Sun but you can do this for any object for which you have the R.A. and Dec.</p>
<p><a href="https://i.sstatic.net/5DvOZ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/5DvOZ.png" alt="Sunlit Earth" /></a></p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt

# for right now, from Google:
RA = (3 + 29/60 + 38/3600) * 360 / 24 # deg for the Sun
dec = +19. # deg for the Sun
GST = (16 + 15/60) * 360 / 24 # deg 

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

latsub = np.radians(dec)
lonsub = np.radians(RA - GST) # where these have already been converted to degrees

theta = np.linspace(0, twopi, 10001)
x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]

tilt = halfpi - np.radians(dec)
x, z = [x * f(tilt) for f in (np.cos, np.sin)]

loncirc = np.arctan2(y, x) - lonsub
loncirc = np.mod(loncirc + pi, twopi) - pi
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))

dloncirc = np.abs(loncirc[1:] - loncirc[:-1])
jump = np.argmax(dloncirc)

latcirc, loncirc = [thing[:-1] for thing in (latcirc, loncirc)]
loncirc[jump] = np.nan
latcirc[jump] = np.nan

locd, lacd = [np.degrees(thing) for thing in (loncirc, latcirc)]
losd, lasd = [np.degrees(thing) for thing in (lonsub, latsub)]
losd = np.mod(losd + 180, 360) - 180 # put it within the maps -180 to +180 range

plt.plot(locd, lacd)
plt.xlim(-182, 182)
plt.ylim(-90, 90)
x, y = 121.5, 25 # Taipei
plt.plot([x], [y], 'or')
plt.text(x, y, ' taipei')
plt.plot([-losd], [lasd], 'oy')
plt.text(-losd, lasd, ' sun')

plt.title('Sunlit Earth 08:40 local time TPE, sunrise was 5:09')
plt.gca().set_aspect('equal')
plt.show()
</code></pre>

CC BY-SA 4.0
2021-05-16T00:30:10.570Z
[I'm 8 hours from Greenwich, why is my sidereal time 10 minutes different from theirs?](https://astronomy.stackexchange.com/q/43819/7982)
uhoh
CC BY-SA 4.0
2021-07-02T11:11:26.077Z
I don't understand this pattern completely.
It's been demonstrated on the Sun. How about the Sirius star for example?
Geographos
CC BY-SA 4.0
2021-07-02T11:17:09.077Z
OK, I got it a bit.
The biggest uncertainty is the declination.
dec = +23. # deg for the Sun    
how will I know that declination is +23,37 ie? It looks like only integers come here. I need the declination more detailed. Could it be possible?
Geographos
CC BY-SA 4.0
2021-07-02T11:19:31.760Z
and how can I define the date and hour here?
Geographos
CC BY-SA 4.0
2021-07-02T12:04:31.820Z
Hi @MKR I see your messages, I will take a look and this again (it's been a while) and then see if/how I can respond...
uhoh
CC BY-SA 4.0
2021-05-16T00:10:07.660Z
2021-05-16T02:33:34.593Z
2021-05-16T02:33:34.593Z
uhoh
<p>I am looking for an option to present the given star visibility world map likewise day &amp; night map or moonlight world map like below:</p>
<p><a href="https://i.sstatic.net/r0j8U.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/r0j8U.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/1mlE7.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1mlE7.png" alt="enter image description here" /></a></p>
<hr />
<p><strong>Question:</strong> Is it possible to generate a map like this for some particular star or at least the point in the sky defined by celestial coordinates?</p>

<night-sky><sky-survey><star-maps><topographic-map>
Star visibility above horizon map likewise day & light map or moonlight world map
Geographos
# Star visibility above horizon map likewise day & light map or moonlight world map

(Asked by: username_0 on 2021-05-15)

I am looking for an option to present the given star visibility world map likewise day & night map or moonlight world map like below:

[![enter image description here](https://i.sstatic.net/r0j8U.png)](https://i.sstatic.net/r0j8U.png)

[![enter image description here](https://i.sstatic.net/1mlE7.png)](https://i.sstatic.net/1mlE7.png)

---

**Question:** Is it possible to generate a map like this for some particular star or at least the point in the sky defined by celestial coordinates?

**Tags:** night-sky, sky-survey, star-maps, topographic-map


---

## 1 answers

### Answer by username_1 on 2021-05-16. Score: 4

As long as we can live with a few tens of kilometers accuracy, we can treat the Earth's surface as a sphere.
The subpoint of the object (the point on Earth at which the object appears at the zenith, overhead, or the central point on Earth seen from the view of the object if they could see Earth clearly) is simply given by
$$\text{lat} = \text{Dec}$$
$$\text{lon} = \text{R.A.} - \text{GST}$$
where GST is Greenwich Sidereal time
To make your plot we draw a great circle around the subpoint, then make a simple Equirectangular projection or lat/lon projection.
In Python for example, calculate the latitude and longitude of the subpoint,
latsub = np.radians(dec)
lonsub = np.radians(RA - GST) # where these have already been converted to degrees

Pick a bunch of points to make a circle
theta = np.linspace(0, 2*np.pi, 10001)

Draw a unit circle around Earth's equator, assuming the subpoint is at the north pole (dec = +90 degrees).
x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]

Tilt it by 90 - dec
tilt = np.pi/2 - np.radians(dec)
x, z = [x * f(tilt) for f in (np.cos, np.sin)]

calculate the latitude and longitude of the tilted circle by shifting the longitude (subtracting the longitude of the subpoint)
loncirc = np.arctan2(y, x) - lonsub
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))

Put it all together and convert units properly, and you get something like this. I did it for the Sun but you can do this for any object for which you have the R.A. and Dec.

import numpy as np
import matplotlib.pyplot as plt

# for right now, from Google:
RA = (3 + 29/60 + 38/3600) * 360 / 24 # deg for the Sun
dec = +19. # deg for the Sun
GST = (16 + 15/60) * 360 / 24 # deg 

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

latsub = np.radians(dec)
lonsub = np.radians(RA - GST) # where these have already been converted to degrees

theta = np.linspace(0, twopi, 10001)
x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]

tilt = halfpi - np.radians(dec)
x, z = [x * f(tilt) for f in (np.cos, np.sin)]

loncirc = np.arctan2(y, x) - lonsub
loncirc = np.mod(loncirc + pi, twopi) - pi
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))

dloncirc = np.abs(loncirc[1:] - loncirc[:-1])
jump = np.argmax(dloncirc)

latcirc, loncirc = [thing[:-1] for thing in (latcirc, loncirc)]
loncirc[jump] = np.nan
latcirc[jump] = np.nan

locd, lacd = [np.degrees(thing) for thing in (loncirc, latcirc)]
losd, lasd = [np.degrees(thing) for thing in (lonsub, latsub)]
losd = np.mod(losd + 180, 360) - 180 # put it within the maps -180 to +180 range

plt.plot(locd, lacd)
plt.xlim(-182, 182)
plt.ylim(-90, 90)
x, y = 121.5, 25 # Taipei
plt.plot([x], [y], 'or')
plt.text(x, y, ' taipei')
plt.plot([-losd], [lasd], 'oy')
plt.text(-losd, lasd, ' sun')

plt.title('Sunlit Earth 08:40 local time TPE, sunrise was 5:09')
plt.gca().set_aspect('equal')
plt.show()

> **Comments:**
> **username_1**: [I'm 8 hours from Greenwich, why is my sidereal time 10 minutes different from theirs?](https://astronomy.stackexchange.com/q/43819/7982)  
> **username_0**: I don't understand this pattern completely. It's been demonstrated on the Sun. How about the Sirius star for example?  
> **username_0**: OK, I got it a bit. The biggest uncertainty is the declination. dec = +23. # deg for the Sun     how will I know that declination is +23,37 ie? It looks like only integers come here. I need the declination more detailed. Could it be possible?  
> **username_0**: and how can I define the date and hour here?  
> **username_1**: Hi @MKR I see your messages, I will take a look and this again (it's been a while) and then see if/how I can respond...  

---

