The project is almost finished

ahmedhosssam

Ahmed Hossam

Posted on July 11, 2024

The project is almost finished

So, after I finished the main refactoring and the documentation strings, I started the testing phase.
I started with the smallest functions that don't call any non-tested function, and then went to the bigger functions.

After finishing the main unit tests, I had a problem, the columns in the input data from the hek have different units, for example, for event_coord1 we have some rows with degrees and others with arcsecond, and the result would be a column with one unit, the degrees are converted into arcseconds or vice versa. So me and my mentor decided to:

  1. Merge event_coord1 and event_coord2 and event_coord3 into one SkyCoord object for every row.
if is_coord_prop:
    event_coord_col = []
    for idx in range(len(table['event_coord1'])):
        unit = get_unit(table['event_coordunit'][idx])
        coord1 = table['event_coord1'][idx]
        coord2 = table['event_coord2'][idx]
        coord3 = table['event_coord3'][idx]
        frame = 'helioprojective' if unit==u.arcsec else 'icrs'
        if coord3:
            event_coord = SkyCoord(coord1*unit, coord2*unit, coord3*unit, frame=frame)
        else:
            event_coord = SkyCoord(coord1*unit, coord2*unit, frame=frame)
        event_coord_col.append(event_coord)
    event_coord_col = Column(event_coord_col, name='event_coord')
    del table['event_coord1']
    del table['event_coord2']
    del table['event_coord3']
    table.add_column(event_coord_col)
Enter fullscreen mode Exit fullscreen mode
  1. The other columns that has multiple units would be converted into astropy column first and then added to the table
if not attribute.get("is_chaincode"):
    new_column = Column(new_column, name=table[attribute["name"]], dtype=u.Quantity)
Enter fullscreen mode Exit fullscreen mode

After this, I started a mission of checking the resulting columns and comparing them to the input and also to the expected results. And here I saw some surprises:
There are columns that don't even exist in the HEK Feature/Event Types definitions
:

  • hgc_coord
  • hpc_coord
  • hgs_coord
  • hrc_coord
  • hgc_boundcc
  • hpc_boundcc
  • hgs_boundcc
  • hrc_boundcc

So I had to add them into coord_properties.json file with double-checking the frames.

hgc_coord, hpc_coord, hgs_coord, hrc_coord are points, so they had to be converted into PointSkyRegion.
And they were supported in parse_chaincode:

def parse_chaincode(value, attribute, unit, time):
.
.
.
    if attribute.get("is_point"):
        coordinates = value.replace("POINT(", "").replace(")", "").split()
        coord_list = [float(coordinate) for coordinate in coordinates]
        coord_list[0] *= coord1_unit
        coord_list[1] *= coord2_unit
        if attribute["frame"] == "heliocentric":
            center_sky = SkyCoord(coord_list[0], coord_list[1], [1]* len(coord_list) * u.AU, obstime=time, observer=observer, representation_type="cylindrical", frame=attribute["frame"])
            return PolygonSkyRegion(vertices=center_sky)
        else:
            center_sky = SkyCoord(coord_list[0], coord_list[1], obstime=time, observer=observer, frame=attribute["frame"])
        return PointSkyRegion(center=center_sky)
.
.
.
Enter fullscreen mode Exit fullscreen mode

So far so good.

Also, some tests in test_hek were simplified, instead of checking if the object is instance of PolygonSkyRegion or PointSkyRegion, we can just check if it's instance of SkyRegion which is the base class for all regions defined in celestial coordinates.

Also, another update to parse_chaincode, the observation time were added to the parameter to be added to all the region objects. We used the column event_starttime to specify the obstime of the event. And we also I added the observer as earth, I still don't know if this would be correct for the different conditions and different queries of the hek, but we will see, and also I wrote a comment to highlight the assumption.

Here is the complete implementation of parse_chaincode:

def parse_chaincode(value, attribute, unit, time):
    """
    Parses a string representation of coordinates and convert them into a PolygonSkyRegion object
    using units based on the specified coordinate frame.

    Parameters
    ----------
    value : str
        A polygon defined using vertices in sky coordinates.
    attribute : dict
        An object from "coord_properties.json"
    unit : str
        The unit of the coordinates
    time : `~astropy.time.core.Time`
        An event_starttime row parsed into astropy time.

    Returns
    -------
    `PolygonSkyRegion`
        A polygon defined using vertices in sky coordinates.

    Raises
    ------
    IndexError
        Because ``value`` does not contain the expected '((' and '))' substrings.
    UnitConversionError
        Because the units set by ``coord1_unit`` or ``coord2_unit`` are incompatible with the values being assigned.
    """
    observer = 'earth' # There is an assumption that earth is the observer.
    coord1_unit = u.deg
    coord2_unit = u.deg
    if attribute["frame"] == "helioprojective":
        coord1_unit = u.arcsec
        coord2_unit = u.arcsec
    elif attribute["frame"] == "heliocentric":
        coord1_unit = u.R_sun # Nominal solar radius
    elif attribute["frame"] == "icrs":
        coord1_unit = get_unit(unit)
        coord2_unit = get_unit(unit)

    if attribute.get("is_point"):
        coordinates = value.replace("POINT(", "").replace(")", "").split()
        coord_list = [float(coordinate) for coordinate in coordinates]
        coord_list[0] *= coord1_unit
        coord_list[1] *= coord2_unit
        if attribute["frame"] == "heliocentric":
            center_sky = SkyCoord(coord_list[0], coord_list[1], [1]* len(coord_list) * u.AU, obstime=time, observer=observer, representation_type="cylindrical", frame=attribute["frame"])
            return PolygonSkyRegion(vertices=center_sky)
        else:
            center_sky = SkyCoord(coord_list[0], coord_list[1], obstime=time, observer=observer, frame=attribute["frame"])
        return PointSkyRegion(center=center_sky)
    coordinates_str = value.split('((')[1].split('))')[0]
    coord1_list = [float(coord.split()[0]) for coord in coordinates_str.split(',')] * coord1_unit
    coord2_list = [float(coord.split()[1]) for coord in coordinates_str.split(',')] * coord2_unit
    vertices = {}
    if attribute["frame"] == "heliocentric":
        vertices = SkyCoord(coord1_list, coord2_list, [1]* len(coord1_list) * u.AU, obstime=time, observer=observer, representation_type="cylindrical", frame="heliocentric")
    else:
        vertices = SkyCoord(coord1_list, coord2_list, obstime=time, observer=observer, frame=attribute["frame"])

    return PolygonSkyRegion(vertices=vertices)
Enter fullscreen mode Exit fullscreen mode

One major outcome of this project was to make using hek and acquiring data from it easy, and this happened when I saw some errors from the CI and when I checked it I saw that one example of overplot_hek_polygon.py should be modified due to the interface changing and the returned data types.
This was the initial using of hek in this example:

ch = responses[response_index]
p1 = ch["hpc_boundcc"][9:-2]
p2 = p1.split(',')
p3 = [v.split(" ") for v in p2]
ch_date = parse_time(ch['event_starttime'])
ch_boundary = responses[response_index]["hpc_boundcc"]

##############################################################################
# The coronal hole was detected at different time than the AIA image was
# taken so we need to rotate it to the map observation time.

ch_boundary = SkyCoord(
    [(float(v[0]), float(v[1])) * u.arcsec for v in p3],
    obstime=ch_date, observer="earth",
    frame=frames.Helioprojective)
rotated_ch_boundary = solar_rotate_coordinate(ch_boundary, time=aia_map.date)
Enter fullscreen mode Exit fullscreen mode

And now it became this:

ch = responses[response_index]
ch_boundary = responses[response_index]["hpc_boundcc"].vertices

# The coronal hole was detected at different time than the AIA image was
# taken so we need to rotate it to the map observation time.

rotated_ch_boundary = solar_rotate_coordinate(ch_boundary, time=aia_map.date)
Enter fullscreen mode Exit fullscreen mode

Yeah, and that's it.
The remaining parts are writing the user documentation, and double-checking the returned data from different events and different queries.

💖 💪 🙅 🚩
ahmedhosssam
Ahmed Hossam

Posted on July 11, 2024

Join Our Newsletter. No Spam, Only the good stuff.

Sign up to receive the latest update from our blog.

Related