Bug or misunderstanding of Conjunction.lla_footprint()

[CassandraAuri]

Describe the issue:

I am following the up to date Conjunction example. That example works with a toy satellite but using empheramis data, I am unsure of how to translate the altitude of the satellite to the footprint altitude.

I have initialized the Conjunction class with conj=asilib.Conjunction(asi, (time, (lat_sat,lon_sat,alt_sat)) then using conj.lla_footprint(alt=110) which returns None and messes up map_azel( see error)

In conjunction_py the comments are out of date so I am unsure what the proper return should be and am quite shocked to see its just one Non. instead an array of length alt_sat of Nones.

Reproduce the code example:

from datetime import datetime, timedelta
import string
import matplotlib.pyplot as plt
import matplotlib.dates
import matplotlib.patches
import matplotlib.gridspec
import numpy as np
import pyaurorax.ephemeris as ae
import asilib
import asilib.asi
time_range = (datetime(2021, 3, 18, 8, 10),
              datetime(2021, 3, 18, 8, 30))

area_box_km = (20, 20)
#skymap_dict = asilib.load_skymap(asi_array_code, location_code, time_range[0])
print(asilib.config['ASI_DATA_DIR'])
global platforms
platforms = ["swarma"]
sky_map_values = ['trex_nir', 'rabb']  
rearth = 6378.1370

def emph():
    print(platforms, "test")

    def coordinates_extrapolation(array):
        n = int((time_range[1] - time_range[0]
                 ).total_seconds() / 3)  # 3 second␣
        latvelocitytotal, lonvelocitytotal = [], []
        cadence = int(60/3)
        # every 3 seconds theres a time component
        time = np.array([time_range[0] + timedelta(seconds=i*3)
                         for i in range(n)])
        for i in range(len(platforms)):
            latvelocity = []
            lonvelocity = []
            for j in range(len(array[0][0])-1):
                latvelocity.append(np.linspace(
                    array[0][i][j], array[0][i][j+1], cadence))
                lonvelocity.append(np.linspace(
                    array[1][i][j], array[1][i][j+1], cadence))
            latvelocitytotal.append(np.reshape(latvelocity, -1))
            lonvelocitytotal.append(np.reshape(lonvelocity, -1))
        print("ree")

        def altitude():
            altitude_total = []
            altepop = []
            altswarma = []
            altswarmb = []
            altvelocityepop = []
            altvelocityswarma = []
            altvelocityswarmb = []
            altswarmc = []
            altvelocityswarmc = []
            for i in range(len(platforms)):

                if(platforms[i] == "epop"):
                    for k in range(len(em1.data)):
                        if(em1.data[k].data_source.platform == "epop"):
                            altepop.append(
                                em1.data[k].metadata["radial_distance"]-rearth)
                    for j in range(len(altepop)-1):
                        altvelocityepop.append(np.linspace(
                            altepop[j], altepop[j+1], 20))
                    altvelocityepop = np.reshape(altvelocityepop, -1)
                    altitude_total.append(altvelocityepop)

                elif(platforms[i] == "swarma"):
                    for k in range(len(em1.data)):
                        if(em1.data[k].data_source.platform == "swarma"):
                            altswarma.append(
                                em1.data[k].metadata["radial_distance"]-rearth)
                    for j in range(len(altswarma)-1):
                        altvelocityswarma.append(np.linspace(
                            altswarma[j], altswarma[j+1], 20))
                    altvelocityswarma = np.reshape(altvelocityswarma, -1)
                    altitude_total.append(altvelocityswarma)

                elif(platforms[i] == "swarmb"):
                    for k in range(len(em1.data)):
                        if(em1.data[k].data_source.platform == "swarmb"):
                            altswarmb.append(
                                em1.data[k].metadata["radial_distance"]-rearth)
                    for j in range(len(altswarmb)-1):
                        altvelocityswarmb.append(np.linspace(
                            altswarmb[j], altswarmb[j+1], 20))
                    altvelocityswarmb = np.reshape(altvelocityswarmb, -1)
                    altitude_total.append(altvelocityswarmb)

                elif(platforms[i] == 'swarmc'):
                    for k in range(len(em1.data)):
                        if(em1.data[k].data_source.platform == "swarmc"):
                            altswarmc.append(
                                em1.data[k].metadata["radial_distance"]-rearth)
                    for j in range(len(altswarmc)-1):
                        altvelocityswarmc.append(np.linspace(
                            altswarmc[j], altswarmc[j+1], 20))
                    altvelocityswarmc = np.reshape(altvelocityswarmc, -1)
                    altitude_total.append(altvelocityswarmc)
                else:
                    pass

            return altitude_total

            # for i in range(len(platforms)):
            # for j in range(len(platforms)):
        altitude_array = altitude()
        return [time, latvelocitytotal, lonvelocitytotal, altitude_array]

    # image_bounds = imager_bounds()
    global em1
    em1 = ae.search(time_range[0], time_range[1],
                    platforms=platforms)
    lattotal, lontotal = [], []
    for k in range(len(platforms)):
        latstarts, lonstarts = [], [],
        for i in range(len(em1.data)):  # 3 spacecraft
            # sees if the data corresponds to the correct space-craft
            if(em1.data[i].data_source.platform == platforms[k]):
                latstarts.append(em1.data[i].location_geo.lat)  # starting
                lonstarts.append((em1.data[i].location_geo.lon))  # ending
        lattotal.append(latstarts)
        lontotal.append(lonstarts)
    return coordinates_extrapolation(
        np.array([lattotal, lontotal]))

 # [time, satellite_lattitute,.] #add everything in main

def main():
    print(__name__)

    save_file = []
    asi_array_code = sky_map_values[0]
    location_code = sky_map_values[1]
    fig, ax = plt.subplots(
        3, 1, figsize=(7, 10), gridspec_kw={'height_ratios': [4, 1, 1]}, constrained_layout=True
    )

    # retrieves emphermaris data in the format [time, [sattelite_longitudes], [satellite_lattitudes], [satelliite_alitiudes]]
    data = emph()
    sat_time = data[0]

    # Finds the footprint of selected region with each selected spacecraft

    def animator():

        for i, (time, image, _, im) in enumerate(movie_generator):
            # Plot the entire satellite track, its current location, and a 20x20 km box
            # around its location.
            ax[1].clear()
            ax[2].clear()
            for j in range(len(sat_azel_pixels_total)):
                ax[0].plot(sat_azel_pixels_total[j][:, 0],
                           sat_azel_pixels_total[j][:, 1], 'blue')
                ax[0].scatter(sat_azel_pixels_total[j][i, 0], sat_azel_pixels_total[j][i, 1],
                              c='red', marker='o', s=50)
                ax[0].contour(area_mask_total[j][i, :, :],
                              levels=[0.99], colors=['yellow'])

                ax[1].plot(sat_time, nearest_pixel_intensity_total[j])
                ax[2].plot(sat_time, area_intensity_total[j])

                # Plot the ASI intensity along the satellite path
            vline1 = ax[1].axvline(time, c='b')
            vline2 = ax[2].axvline(time, c='b')

            # Annotate the location_code and satellite info in the top-left corner.
            location_code_str = (
                f'{asi_array_code}/{location_code} '
                f'LLA=({asi.meta["lat"]:.2f}, '
                f'{asi.meta["lon"]:.2f}, {asi.meta["alt"]:.2f})'
            )
            text_obj = ax[0].text(
                0,
                1,
                location_code_str,
                va='top',
                transform=ax[0].transAxes,
                color='red',
            )
            ax[1].set(ylabel='ASI intensity\nnearest pixel [counts]')
            ax[2].set(xlabel='Time',
                      ylabel='ASI intensity\n10x10 km area [counts]')

    def ASI_logic():

        alt = 110  # footprint value
        if(asi_array_code.lower() == 'themis'):
            asi = asilib.asi.themis(
                location_code, time_range=time_range, alt=alt)
        elif(asi_array_code.lower() == 'rego'):
            asi = asilib.asi.rego(
                location_code, time_range=time_range, alt=alt)
            print("test")
        elif(asi_array_code.lower() == 'trex_nir'):
            asi = asilib.asi.trex.trex_nir(
                location_code, time_range=time_range, alt=alt)
        return asi
    asi = ASI_logic()
    # Initiate the movie generator function. Any errors with the data will be␣

    movie_generator = asi.animate_fisheye_gen(  # initaliziation
        ax=ax[0], azel_contours=True, overwrite=True, cardinal_directions='NE'
    )
    # Use the generator to get the images and time stamps to estimate mean the ASI
    # brightness along the satellite path and in a (20x20 km) box.
    sat_azel_pixels_total, nearest_pixel_intensity_total, area_intensity_total, area_mask_total = [], [
    ], [], []  # Creates empty arrays for data from satellites to go, which is there rendered in animator

    def Animation_logic_generator():
        nonlocal sat_azel_pixels_total, nearest_pixel_intensity_total, area_intensity_total, area_mask_total

        for i in range(len(platforms)):  # length of spacecraft

            conjunction_obj = asilib.Conjunction(asi, (sat_time, np.array(  # lat_sattelite, long, alt
                [data[1][i], data[2][i], data[3][i]]).T))
            print(data[3][i])
            blah = conjunction_obj.lla_footprint(alt=110)
            print(blah)
            # Normally the satellite time stamps are not the same as the ASI.
            # You may need to call Conjunction.interp_sat() to find the LLA coordinates
            # at the ASI timestamps.
            # Map the satellite track to the imager's azimuth and elevation coordinates and
            # image pixels. NOTE: the mapping is not along the magnetic field lines! You need
            # to install IRBEM and then use conjunction.lla_footprint() before
            # calling conjunction_obj.map_azel.
            sat_azel, sat_azel_pixels = conjunction_obj.map_azel()
            nearest_pixel_intensity = conjunction_obj.intensity(box=None)
            area_intensity = conjunction_obj.intensity(box=(10, 10))
            area_mask = conjunction_obj.equal_area(box=(10, 10))

            # Need to change masked NaNs to 0s so we can plot the rectangular area contours.
            area_mask[np.where(np.isnan(area_mask))] = 0
            sat_azel_pixels_total.append(sat_azel_pixels)
            nearest_pixel_intensity_total.append(nearest_pixel_intensity)
            area_intensity_total.append(area_intensity)
            area_mask_total.append(area_mask)
        # sat_azel_pixels, area_box_mask_2, asi_brightness_2 for each satellite
    Animation_logic_generator()
    animator()

    print(
        f'Movie saved in {asilib.config["ASI_DATA_DIR"] / "animations"}', )

    movie_container = 'mp4'
    movie_address = f'{time_range[0].strftime("%Y%m%d_%H%M%S")}_' \
        f'{time_range[1].strftime("%H%M%S")}_' \
        f'{asi_array_code.lower()}_{location_code.lower()}_fisheye.{movie_container}'

    movie_address_total = asilib.config["ASI_DATA_DIR"] / \
        'animations'/movie_address
    print(movie_address_total)
    save_file.append(movie_address_total)
    return save_file

if __name__ == '__main__':  # arg parse look up
    test = main()

Error message:

Traceback (most recent call last):
  File "c:/Users/1101w/Clone/Programming-8/Summer 2023/AnimationSattelite.py", line 262, in <module>
    test = main()
  File "c:/Users/1101w/Clone/Programming-8/Summer 2023/AnimationSattelite.py", line 243, in main
    Animation_logic_generator()
  File "c:/Users/1101w/Clone/Programming-8/Summer 2023/AnimationSattelite.py", line 231, in Animation_logic_generator
    sat_azel, sat_azel_pixels = conjunction_obj.map_azel()
  File "C:\Users\1101w\anaconda3\lib\site-packages\asilib\conjunction.py", line 315, in map_azel
    az, el, _ = pymap3d.geodetic2aer(
  File "C:\Users\1101w\anaconda3\lib\site-packages\pymap3d\aer.py", line 110, in geodetic2aer
    e, n, u = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
  File "C:\Users\1101w\anaconda3\lib\site-packages\pymap3d\enu.py", line 199, in geodetic2enu
    x1, y1, z1 = geodetic2ecef(lat, lon, h, ell, deg=deg)
  File "C:\Users\1101w\anaconda3\lib\site-packages\pymap3d\ecef.py", line 68, in geodetic2ecef
    lat, ell = sanitize(lat, ell, deg)
  File "C:\Users\1101w\anaconda3\lib\site-packages\pymap3d\utils.py", line 62, in sanitize
    raise ValueError("-pi/2 <= latitude <= pi/2")
ValueError: -pi/2 <= latitude <= pi/2
PS C:\Users\1101w\Clone\Programming-8>

Runtime information:

0.16.0 3.8.3 (default, Jul 2 2020, 17:30:36) [MSC v.1916 64 bit (AMD64)] Windows-10-10.0.19041-SP0

Context for the issue:

No response

[mshumko]

My apologies for the late reply @CassandraAuri. I'm traveling these next two weeks and I haven't had the chance to run your code yet.

The conjunction_obj.lla_footprint method does not return anything. It does, however, modify the internal conjunction_obj.sat pandas dataframe attribute.

The exception is raised my pymap3d which does the actual mapping. My first though is that there is an improper latitude value(s) passed into it. Can you verify that a) max(lat) < 90 deg, min(lat) > -90 deg, and there are no NaN or -1E31 values? The -1E31 is IRBEM's invalid value flag. You will need to check the values in conjunction_obj.sat right before you call map_azel().

[CassandraAuri]

So my altitude array was just empty as I was revamping it so the program obviously wouldn't know where my satellite was, otherwise it seems to be working as intended and I get the exact same image as pre-imager class implementation.

Cassandra