Limite jour/nuit sur une carte

À la demande de Thierry Loiseau (astrophoto.free.fr, merci pour le défi lancé !), j’ai rédigé un script Python qui trace la limite jour/nuit sur une carte du monde, pour une date donnée. Le point jaune indique la position survolée par le Soleil.

Évolution au cours d’une journée :

Évolution sur une année :

Le script Python :

Pour ce script, j’utilise les modules habituels (matplotlib, numpy, datetime) ainsi que le module basemap qui permet de tracer des cartes dans un graphique, et que j’ai découvert pour l’occasion.

Module basemap

Avant de choisir l’option basemap, j’ai suivi la piste d’insérer une carte en image de fond sur le graphique. Je savais qu’il ne serait pas simple de placer les coordonnées, même avec une carte de type Mercator. J’ai trouvé sur le web une expression littérale de la latitude d’un point sur la carte, en fonction de sa latitude sur terre, mais cela ne collait pas du tout.

C’est alors que je suis tombé sur le module basemap, qui crée un « axe » (au sens de matplotlib) avec une carte assez facile à personnaliser. Placer sur cette carte les coordonnées calculées (position du Soleil, points calculés de la courbe) n’allaient pas de soi, mais en tâtonnant on trouve assez vite comment utiliser ce module.

Créer et personnaliser la carte est facile, avec les exemples de scripts fournis dans la documentation.

On peut choisir le type de projection, la couleur de remplissage des océans et des continents, frontières… L’option « bluemarble » que j’ai utilisée dessine un fond de carte de type image satellite. La résolution de cette image se choisit avec le paramètre « scale » ; attention pour une forte résolution l’exécution du script est longue, et l’image finale assez lourde.

#  création de la carte :
m = Basemap(
    projection='mill', # type de projection (voir documentation Basemap)
    llcrnrlon=-180,urcrnrlon=180,lat_ts=0, resolution='c',
    llcrnrlat=-75, urcrnrlat=75,  # limites en latitude
            )
m.drawcountries(linewidth = 0.3, color='white', zorder=1)
m.bluemarble(scale = 0.5)

On peut choisir une option plus légère :

#  création de la carte :
m = Basemap(
    projection='mill', # type de projection (voir documentation Basemap)
    llcrnrlon=-180,urcrnrlon=180,lat_ts=0, resolution='c',
    llcrnrlat=-75, urcrnrlat=75,  # limites en latitude
            )
m.drawcountries(linewidth = 0.3, color='white', zorder=1)
m.drawcoastlines(linewidth = 0.4)
m.fillcontinents(color='peru')
#  ------------------------------------------------
# tracé des méridiens et parallèles :
m.drawmeridians(np.arange(-180., 180., 15),linewidth=0.2,color='gray',
                zorder = 1,
                dashes = [1, 0],  # trait continu
                )
m.drawparallels(np.arange(-75, 75, 15),linewidth=0.2,color='gray',
                dashes=[1,0],  # trait continu
                zorder=1,
                )

Il est très facile de choisir de tracer les méridiens et parallèles, de choisir leur espacement et leur mise en forme.

En revanche, une fois la carte personnalisée, quand j’ai voulu y superposer les courbes et position du Soleil, ça a été moins évident… Pour que l’axe de la carte soit utilisable pour les commandes usuelles de matplotlib, il faut le créer de cette manière :

fig = plt.figure(tight_layout=True)
ax = fig.add_subplot(111)

Cela ne fonctionne pas avec :

fig = plt.figure(tight_layout=True)
ax = plt.subplot(111)

Ensuite, il faut comprendre que les coordonnées sur la carte doivent être calculées par une commande du module basemap. Si on fait ax.scatter(0,50), on ne placera pas un point à la longitude 0 et la latitude 50. La méthode est la suivante :

x = longitude
y = lat_deg
x, y = m(x, y)  # transformation des coordonnées
ax.plot(x, y)  # tracé sur la carte

Il faut faire de même si on veut tracer les tropiques, ou si on veut la coordonnée correspondant à la latitude 90°, par exemple.

Calculs astronomiques

Je pensais les calculs compliqués au départ. Voici la procédure :

• on crée une liste de valeur de longitudes, de -180 à 180°.
• on définit l’angle horaire du Soleil au méridien de Greenwich (d’après la variable d’heure définie au début du script).
• pour chaque longitude, on calcule l’angle horaire du Soleil qui lui correspond.
• enfin, pour chaque longitude, on calcule la latitude pour laquelle le Soleil est sur l’horizon à cet angle horaire. C’est plus simple qu’on ne le pense. La formule donnant la hauteur du Soleil h en fonction de sa déclinaison d, de son angle horaire H et de sa latitude L s’écrit sin(h) = sin(d) * sin(L) + cos(d)*cos(L)*cos(H). On veut h=0 (Soleil sur l’horizon, hauteur nulle), H et d sont déjà définis. Il suffit de modifier la formule pour exprimer la latitude L : cos(H) = – tan(d)*tan(L) (on retrouve la formule qui donne la durée du jour, au passage). Donc tan(L) = -cos(H) / tan(d)