Ground track using cartopy

[gregoriomarchesini]

Hi I have developed a small piece of code that can visualize ground tracks using cartopy The available list of map is listed in the documentation of the function.

I tried it and it works for me. You have to install cartopy from conda, but nothing else is required (apart from numpy, but this is pretty normal to have).

If you think it is useful, let me know. I hope it will help

Have a super day!

def groundtrack_generator(settings    :dict = None,
                          groundtracks :list = []) :

    """
    Parameters 
    ----------

    settings   (dict)     : define settings for the ground track plot (see intialisation
                           of settings for a complete list of settings). If an empty 
                           dictionary is given as input, the defoult settigns are used 

    groundtracks (list)  : define list of ground tracks to plot. Each ground track must be 
                           defined as a dictionary :

                           groundtrack  = {
                               "name" :
                               "lonlat" : np.ndarray[[N,2]] defining the lat and long of the 
                                          ground track. col[0] = longitude col[1] = latitude.
                                          Units are radiants.  
                                          }

    Returns
    --------
    only figure of the ground track is returned

    Description :
    -------------
    plots a set of ground tracks over a specified map. 

    have a look at :
    https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html

    for all maps projections and types

    Maps

    1  PlateCarree
    2  Mollweide
    3  AlbersEqualArea
    4  AzimuthalEquidistant
    5  LambertConformal
    6  LambertCylindrical
    7  Mercator
    8  Miller
    9  Mollweide
    10 Orthographic
    11 Robinson
    12 Sinusoidal
    13 Stereographic
    14 TransverseMercator
    15 InterruptedGoodeHomolosine
    16 RotatedPole

    ## settings

    _settings          = {"central_long" : 0,                                 # central longitude
                          "color"        : [(0.0,0.0,0.0)]*len(groundtracks), # color per ground track
                          "projection"   : "PlateCarree",                     # projection (see list above)
                          "linewidth"    : 1,                                 # line width
                          "grid"         : False,      # grid on the map
                          "grid_linewidth":  1,
                          "grid_color"    : "gray",
                          "grid_alpha"    : 0.3}

    """

    rad2deg = 180/np.pi

    _settings          = {"central_long" : 0,      
                          "color"        : [(0.0,0.0,0.0)]*len(groundtracks),
                          "projection"   : "PlateCarree",
                          "linewidth"    : 1,
                          "grid"         : False,
                          "grid_linewidth":  1,
                          "grid_color"    : "gray",
                          "grid_alpha"    : 0.3}

    if not settings == None :
        for key in settings.keys():
            if key not in _settings.keys() :
                raise ValueError("unknown argument {}".format(key))

            # impose desired settings

    _settings.update(settings)

    try :
        map_projection = getattr(cartopy.crs,_settings["projection"]) 
    except  AttributeError() :
        raise ValueError("unknown map type {}.\n Check again the correct name of the map you want to apply".format(key))

    # start defining the figure
    figure        = plt.figure(figsize=(8,6))
    ax            = plt.axes(projection=map_projection())

    for n,groundtrack in enumerate(groundtracks) :

        plt.scatter(groundtrack["lonlat"][:,0]*rad2deg,groundtrack["lonlat"][:,1]*rad2deg,
                    color     =  _settings["color"][n],
                    s = _settings["linewidth"],
                    label     = groundtrack["name"],
                    transform=cartopy.crs.Geodetic())
        plt.scatter(groundtrack["lonlat"][0,0]*rad2deg,groundtrack["lonlat"][0,1]*rad2deg,
                    s =_settings["linewidth"]*4,
                    color     =  _settings["color"][n],
                    transform=cartopy.crs.Geodetic())

    ax.set_global()

    if _settings["grid"] :
       ax.gridlines(color     =_settings["grid_color"],
                    linestyle ="--",
                    draw_labels=True,
                    linewidth =_settings["grid_linewidth"],
                    alpha     =_settings["grid_alpha"])

    ax.coastlines()
    try :
        ax.stock_img()
    except :
        pass

    ax.legend(bbox_to_anchor =(0.0,0.0))

    plt.show()

[DominicDirkx]

Hi Greg, we had a long weekend holiday here, hence the late reply. This looks like a nice piece of functionality again :) One question on the 'settings' input. Wouldn't it be easier to have each of the settings entries be an explicit input (with default) for the function itself? That way:

• Autocomplete will show which options there are
• Autocomplete will show the default values
• Documenting this on our API will be easier

But, as you know, I am not that good a Python programmer, nor am I very familiar with common practices. Let me know what you think about the above suggestion

[gregoriomarchesini]

Absolutely

I will correct it soon! It is indeed a good suggestion

Super I will repost here the updated version

[gregoriomarchesini]

Hi again

I have corrected the function and I believe now it could be more useful. I inserted the parameters as you asked in the function signature

I have also created a minimal working example so that you can test if everything is working as you wish

import numpy as np
from matplotlib import pyplot as plt

from tudatpy.io import save2txt
from tudatpy.util import result2array
from tudatpy.kernel import constants
from tudatpy.kernel.interface import spice
from tudatpy.kernel import numerical_simulation # KC: newly added by me
from tudatpy.kernel.numerical_simulation import environment_setup
from tudatpy.kernel.numerical_simulation import propagation_setup
from scipy.spatial.transform import Rotation
import cartopy

def split_history(state_history:dict,n_bodies:int):

    """
    Parameters
    -----------

    state_history        (dict ) :  dictionary containing the state_history
                                    of N propagated bodies.
                                    The dictionary is defined as a standard 
                                    state_history following the TUdat convention. 

                                    the keys of the dictionay are the time stamps 
                                    corresponding to each state and the argument is 
                                    the state of the body at the given time stamp

    n_rotation         (int)     : number of propagated bodies

    Returns
    -----------

    state_hostory_book (list)   : list of state_history dictionaries for each propagated body.
                                  The order of the list is given directly by the order 
                                  in which the bodies are propagated.

    Descriptiom
    -----------
    Creates a list of state histories based on the unified state_history
    from from the propagation of N satellites

    """

    state_history_book = []
    time        = [key for key in state_history.keys() ]
    n_variables        = len(state_history[time[0]])
    step               = int(n_variables/n_bodies)

    if  n_variables%n_bodies :
        print(" ERROR : Not possible to split. Number of variables not divisible by number of bodies")
        raise 

    for n in range(n_bodies) :
        current_history = dict()

        for t in state_history.keys():
           current_history[t]= state_history[t][step*n:step*n+step]

        state_history_book.append(current_history)

    return state_history_book

def groundtrack_generator(groundtracks_list :list  = [],
                          labels            :list = [],
                          projection        :str   = "PlateCarree",
                          colors            :list  = [],
                          linewidth         :float = 1 ,
                          grid              :bool  = True,
                          grid_linewidth    :float = 1,
                          grid_color        :str   = "gray",
                          grid_alpha        :float =   0.3) :

    """
    Parameters 
    ----------

    groundtracks_list (List)  : List of numpy.ndarray[[N,2]] defining the latitude and longitude of each groundtrack
                                col[0] = longitude 
                                col[1] = latitude.
                                Units are radiants.  

    labels            (List ) : label per groundtrack
    projection        (Str  ) : projection to be used (see available projections below)
    colors            (List ) : list of colours per groundtrack
    linewidth         (Float) : linewidth of the groundtracks 
    grid              (Bool ) : activates grid on map, defining latitude and longitude reference
    grid_linewidth    (Float) : grid width
    grid_color        (Str  ) : grid color (refer to https://matplotlib.org/3.5.0/tutorials/colors/colors.html#:~:text=Single%20character%20shorthand%20notation%20for%20some%20basic%20colors. for a list of colors)
    grid_alpha        (Float) : grid alpha

    Returns
    --------
    Returns plot of a given set of ground tracks as definied by the groundtrack_list

    Description :
    -------------
    plots a set of ground tracks over a specified map. 

    have a look at :
    https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html
    for all maps projections and types

    Maps

    1  PlateCarree
    2  Mollweide
    3  AlbersEqualArea
    4  AzimuthalEquidistant
    5  LambertConformal
    6  LambertCylindrical
    7  Mercator
    8  Miller
    9  Mollweide
    10 Orthographic
    11 Robinson
    12 Sinusoidal
    13 Stereographic
    14 TransverseMercator
    15 InterruptedGoodeHomolosine
    16 RotatedPole

    """
    try :
        map_projection = getattr(cartopy.crs,projection) 
    except  AttributeError :
        raise ValueError("unknown map type {}.\n Check again the correct name of the map you want to apply".format(projection))

    # start defining the figure
    figure        = plt.figure(figsize=(8,6))
    ax            = plt.axes(projection=map_projection())

    rad2deg = 180/np.pi
    if len(colors) == 0 :
        colors = [(.0,.0,.0,1.)]*len(groundtracks_list)
    elif len(colors) != len(groundtracks_list):
        raise ValueError("colors and groundtracks_list have different lenghts")
    if len(labels) == 0 :
        labels = ["sat{}".format(n) for n in range(len(groundtracks_list))]
    elif len(labels) != len(groundtracks_list):
        raise ValueError("labels sand groundtracks_list have different lenghts")

    for groundtrack,color,label in zip(groundtracks_list,colors,labels) :

        # plot start position with bigger width

        long = groundtrack[:,0]*rad2deg
        lat  = groundtrack[:,1]*rad2deg

        plt.plot(long,lat,
                    color     =  color,
                    linewidth = linewidth,
                    label     = label,
                    transform=cartopy.crs.Geodetic())

        plt.scatter(long[0],lat[0],
                    linewidth=linewidth*5,
                    color     =  color,
                    transform=cartopy.crs.Geodetic(),
                    marker = "*")

    ax.set_global()

    if grid :
       ax.gridlines(color     =grid_color,
                    linestyle ="--",
                    draw_labels=True,
                    linewidth =grid_linewidth,
                    alpha     =grid_alpha)

    try :
        ax.stock_img()
    except :
        pass

    ax.legend()

# define useful constants

julian_day = constants.JULIAN_DAY
hours      = 60*60
minutes    = 60 
rad2deg    = 180/np.pi

# start simulation day 20/10/1998
# simulate for 10 days
simulation_start_epoch = 2451106.50000
simulation_end_epoch   = simulation_start_epoch + 4*hours

# Load SPICE.
spice.load_standard_kernels() 

# Create settings for celestial bodies
bodies_to_create         = ['Earth']      # this must have a list of all the planets to create
global_frame_origin      = 'Earth'        # this is the origin of the refernce system
global_frame_orientation = 'J2000'   # orinetation of the reference system
body_settings = environment_setup.get_default_body_settings(
    bodies_to_create, global_frame_origin, global_frame_orientation) # body settings taken from SPICE.

# Create environment
bodies = environment_setup.create_system_of_bodies(body_settings)

# define constellation of satellites
number_of_satellites = 4
satellites_names     = ["sat{}".format(str(num)) for num in range(number_of_satellites)]

for name in satellites_names : 
    bodies.create_empty_body(name)
    # Set mass of vehicle
    bodies.get_body(name).mass = 500.0

###########################################################################
# CREATE ACCELERATIONS ####################################################
###########################################################################

# Define bodies that are propagated, and their central bodies of propagation.
bodies_to_propagate = satellites_names
central_bodies      = ['Earth']*len(bodies_to_propagate)   # bodies which are contributing with the full gravity and are not only perturbation

# Define accelerations acting on vehicle.
acceleration_settings_satellite = dict(
    Earth=[propagation_setup.acceleration.spherical_harmonic_gravity(2,2)]  
)

global_acceleration_settings = {name: acceleration_settings_satellite for name in satellites_names}

# Create acceleration models.
acceleration_models = propagation_setup.create_acceleration_models(
        bodies, global_acceleration_settings, bodies_to_propagate, central_bodies)

# Define initial state.
standard_initial_state =  np.array([5.194681552632270E+06,-1.833738414039049E+06,4.172975147134797E+06,
                                   4.556144583672027E+03,4.992919206236066E+03,-3.471626733390986E+03])[:,np.newaxis]

phi_intervals            = np.linspace(0,1.9*np.pi,number_of_satellites)
multidimenasional_state  = []
dependent_variables      = []

for phi,name in zip(phi_intervals,satellites_names) :
    rot_mat        = Rotation.from_rotvec(np.array([0,0,1])*phi).as_matrix()
    full_state_rot = np.eye(6)
    full_state_rot[:3,:3] = rot_mat
    full_state_rot[3:,3:] = rot_mat
    updated_state  = (full_state_rot @standard_initial_state)
    multidimenasional_state.append(updated_state)

    dependent_variables.append(propagation_setup.dependent_variable.longitude(name,'Earth'))
    dependent_variables.append(propagation_setup.dependent_variable.geodetic_latitude(name,'Earth'))

multidimenasional_state = np.concatenate(tuple(multidimenasional_state),axis=0)

# Create propagation settings.
termination_settings = propagation_setup.propagator.time_termination( simulation_end_epoch )
propagator_settings  = propagation_setup.propagator.translational(
    central_bodies,
    acceleration_models,
    bodies_to_propagate,
    multidimenasional_state,
    termination_settings,
    output_variables = dependent_variables
)

# Create numerical integrator settings.
fixed_step_size = 10.0
integrator_settings = propagation_setup.integrator.runge_kutta_4(
    simulation_start_epoch,
    fixed_step_size
)

# Create simulation object and propagate dynamics.
dynamics_simulator = numerical_simulation.SingleArcSimulator(
    bodies, integrator_settings, propagator_settings )

# create list of state dictionaries for each satellite

state_history_book            = split_history(dynamics_simulator.state_history,number_of_satellites) #list of state histories
state_history_array_book      = [result2array(state_history) for state_history in state_history_book]
dependat_variables_array      = result2array(dynamics_simulator.dependent_variable_history)

# divide long/lat info from rotation matrix
longlat_coordinates     = dependat_variables_array[::20,1:number_of_satellites*2+1]
groundtracks_list       = [longlat_coordinates[:,ii:ii+2] for ii in range(0,number_of_satellites*2,2)]
labels                  = ["sat{}".format(n) for n in range(len(groundtracks_list))]
color_list              = [(.0,.0,.0,1.) for ii in range(number_of_satellites)]
map_settings = {"projection" : "Sinusoidal",
                "color"      : color_list,
                "grid"       : True,
                "grid_alpha" : 0}

test_maps =  ["Sinusoidal",
              "Stereographic",
             "LambertCylindrical",
             "Mercator",
             "Miller",
             "Mollweide",
             "Orthographic",
             "InterruptedGoodeHomolosine"]

for map in test_maps :

    groundtrack_generator(groundtracks_list = groundtracks_list,
                        colors            = color_list,
                        labels            = labels,
                        projection        = map,
                        linewidth         = 1 ,
                        grid              = True,
                        grid_linewidth    = 1,
                        grid_color        = "gray",
                        grid_alpha        =  0.3)

plt.show()

In addition to Tudat, the packages scipy and cartopy are required. But nothing more special than this

I hope the function will be useful in the future, and if not you can just ignore this current issue.

Thank you for your time and attention Have a nice evening

Greg

[DominicDirkx]

Hi Greg, somehow we lost track of this code your wrote, our bad! I think we can split this in two:

• Add your new Python functions into Tudatpy
• Add the functionality into an example application. Perhaps this one: https://github.com/tudat-team/tudatpy-examples/issues/14 would be a good choice

[gregoriomarchesini]

Dear Dominic

I will do my best to find this piece of code and restore it. It has probably been an year since the last time I opened it. This week I will try to conclude it. Unfortunately my week is super busy until Friday. In the weekend I can try to fix it and see if the code is still working.

Thank you for your time and attention Best Greg

On 14 Dec 2022, at 14:20, DominicDirkx @.***> wrote:

Hi Greg, somehow we lost track of this code your wrote, our bad! I think we can split this in two:

Add your new Python functions into Tudatpy Add the functionality into an example application. Perhaps this one: tudat-team/tudatpy-examples#14 https://github.com/tudat-team/tudatpy-examples/issues/14 would be a good choice — Reply to this email directly, view it on GitHub https://github.com/tudat-team/tudat-space/issues/66#issuecomment-1351343719, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQ76ADECCSZWSF5CH55IDYTWNHCTTANCNFSM5TP2ZN7A. You are receiving this because you authored the thread.

[DominicDirkx]

Hi Greg, no need to hurry! I didn't mean to ask you to do this, I think we can merge your new Python functions into Tudat already (if you want 'credit' in the tudatpy commit history, you can open a pull request, though). The use in the Galileo example will come when we properly write that example.

[gregoriomarchesini]

Hey Dominic

Ok perfect. There is no need for the pull request. I am very happy that my code gets used and I hope it will serve its purpose. In case you need help, I will do my best to help you in time. Super thank you and have a great day

Best Greg

On 15 Dec 2022, at 15:37, DominicDirkx @.***> wrote:

Hi Greg, no need to hurry! I didn't mean to ask you to do this, I think we can merge your new Python functions into Tudat already (if you want 'credit' in the tudatpy commit history, you can open a pull request, though). The use in the Galileo example will come when we properly write that example.

— Reply to this email directly, view it on GitHub https://github.com/tudat-team/tudat-space/issues/66#issuecomment-1353193316, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQ76ADCK7SCXEUWGGXDU5YLWNMULBANCNFSM5TP2ZN7A. You are receiving this because you authored the thread.