DEV Community: Ahmed Hossam

The project is almost finished

Ahmed Hossam — Thu, 11 Jul 2024 19:38:21 +0000

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:

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)

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)

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)
.
.
.

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)

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)

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)

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.

Finishing implementation and doc strings

Ahmed Hossam — Mon, 17 Jun 2024 04:44:45 +0000

The main implementation and documentation strings phase has been finished.

One of the functions that have been deleted was parse_unit, it almost wasn't triggered by any event of HEK events, so I decided to delete it, and we will see in the future if it has any major effects on the code.

Also parse_columns_to_table has been refactored by adding more variables to make the code more readable.

This is the new implementation:

def parse_columns_to_table(table, attributes, is_coord_prop = False):
    """
    Parses the columns in an Astropy table and convert the values into Astropy objects.

    Parameters
    ----------
    table: astropy.table
        Astropy table.
    attributes: list
        A list of HEK unit attributes or coordinate attributes.
    is_coord_prop: bool
        To specify if `attributes` is a list of unit attributes or coordinate attributes.

    Returns
    -------
    `astropy.table`

    Raises
    ------
    TypeError
        If `table` is not an Astropy table.
    KeyError
        If any of the attribute dictionaries are missing required keys (i.e. "name", "unit_prop").
    """
    for attribute in attributes:
        if attribute["name"] in table.colnames and ("unit_prop" in attribute or attribute.get("is_chaincode", False)) and attribute.get("is_unit_prop", True):
            unit_attr = ""
            if is_coord_prop:
                unit_attr = "event_coordunit"
            else:
                unit_attr = attribute["unit_prop"]

            new_column = []
            for idx, value in enumerate(table[attribute["name"]]):
                new_value = ""
                if value in ["", None]:
                    new_value = value
                elif attribute.get("is_chaincode", False):
                    new_value = parse_chaincode(value, attribute, table[attribute["unit_prop"]][idx])
                else:
                    unit = get_unit(table[unit_attr][idx])
                    new_value = value * unit
                new_column.append(new_value)

            table[attribute["name"]] = new_column

    for attribute in attributes:
        if attribute.get("is_unit_prop", False) and attribute["name"] in table.colnames:
            del table[attribute["name"]]
    return table

Also a public interface __all__ has been added to util.py, until now, the public "api" for util.py is:

__all__ = [
    'freeze',
    'parse_times',
    'parse_values_to_quantities',
    'UNIT_FILE_PATH',
    'COORD_FILE_PATH'
]

UNIT_FILE_PATH and COORD_FILE_PATH are used in test_hek.py

There was a problem in get_unit because of some units that wasn't parsed correctly, one of them is ergs per cubic centimeter, because of the first implementation of get_unit:

with u.add_enabled_units([cm2, m2, m3]), u.set_enabled_aliases(aliases):
        # Units for coordinate frames have more than one unit, otherwise it will be just one unit.
        # There is an assumption that coord1_unit, coord2_unit and coord3_unit are the same.
        units = re.split(r'[, ]', unit)
        return u.Unit(units[0].lower())

The function was taking the first "word" of the input and returns the unit, obviously this will be wrong with ergs per cubic centimeter because it will take only ergs and returns u.Unit('erg').
But this is the only case with HEK units, the other units works just fine.
This is the new implementation to correct this behavior:

def get_unit(unit):
    """
    Converts string into astropy unit.

    Parameters
    ----------
    unit: str
        The targeted unit

    Returns
    -------
    unit
        Astropy unit object (e.g. <class 'astropy.units.core.Unit'> or <class 'astropy.units.core.CompositeUnit'>)

    Raises
    ------
    ValueError
        Because `unit` did not parse as unit.

    Notes
    ----
    For the complete list of HEK parameters: https://www.lmsal.com/hek/VOEvent_Spec.html

    """
    cm2 = u.def_unit("cm2", u.cm**3)
    m2 = u.def_unit("m2", u.m**2)
    m3 = u.def_unit("m3", u.m**3)
    erg_per_cm3 = u.def_unit("ergs/cm^3", u.erg/u.ml)

    aliases = {
        "steradian": u.sr,
        "arcseconds": u.arcsec,
        "degrees": u.deg,
        "sec": u.s,
        "emx": u.Mx,
        "amperes": u.A,
        "ergs": u.erg,
        "cubic centimeter": u.ml,
        "square centimeter": cm2,
        "cubic meter": m3,
        "square meter": m2,
        "ergs per cubic centimeter": erg_per_cm3,
    }
    with u.add_enabled_units([cm2, m2, m3]), u.set_enabled_aliases(aliases), warnings.catch_warnings():
        # Units for coordinate frames have more than one unit, otherwise it will be just one unit.
        # There is an assumption that coord1_unit, coord2_unit and coord3_unit are the same.
        warnings.filterwarnings("ignore", category=u.UnitsWarning)
        if unit in aliases:
            unit = u.Unit(aliases[unit])
        else:
            unit = u.Unit(re.split(r'[, ]', unit)[0].lower())
        return unit

Also, the units warnings were ignored because it's irrelevant to the user of HEK.

So, the next phase is testing, we will see if our assumptions about the unit parsing and the other refactoring were right or not.

Good start for the project

Ahmed Hossam — Mon, 27 May 2024 13:23:44 +0000

First 2 weeks summarize:

This week I began working on refactoring hek.py functions.
I started by migrating the finished work in GSoC2023 to a new PR to start working on it.
My first contribution was creating util.py file to include all utility functions needed for hek.py, a lot of functions that was added in HEKClient at first didn't make sense to remain there.

Now the new util.py file includes:

def parse_times(table)
def parse_values_to_quantities(table)
def parse_columns_to_table(table, attributes, is_coord_prop = False)
def parse_unit(table, attribute, is_coord_prop = False)
def parse_chaincode(value, attribute, unit)
def get_unit(unit)

get_unit has been simplified in terms of implementation and interface, this was the first version:

    def get_unit(attribute, str):
        if attribute["is_coord_prop"]:
            coord1_unit, coord2_unit, coord3_unit = None, None, None
            coord_units = re.split(r'[, ]', str)
            if len(coord_units) == 1: # deg
               coord1_unit = coord2_unit = u.Unit(coord_units[0])
            elif len(coord_units) == 2:
                coord1_unit = u.Unit(coord_units[0])
                coord2_unit = u.Unit(coord_units[1])
            else:
                coord1_unit = u.Unit(coord_units[0])
                coord2_unit = u.Unit(coord_units[1])
                coord3_unit = u.Unit(coord_units[2])
            return locals()[attribute["unit_prop"]]
        else:
            return u.Unit(str)

The first thing that has been done is to use unit aliases inside the function using context manager instead of putting the aliases globally.

The whole goal of this function is to parse a string into an astropy unit, but the big part of the function was splitting the string into more than one unit if the input was coordinate units, and then returning the unit assigned to unit_prop. I decided to just remove all of this and convert the unit into an array and return the first index, like this:

   units = re.split(r'[, ]', unit)
   return u.Unit(units[0].lower())

And actually it works just fine with all HEK features and events, so I will keep it like this until some strange error appears.
And also the interface has been simplified to just take the string of the targeted unit.

This is the current version of get_unit:

def get_unit(unit):
    """
    Converts string into astropy unit.

    Parameters
    ----------
    unit: str
        The targeted unit

    Returns
    -------
    unit
        Astropy unit object (e.g. <class 'astropy.units.core.Unit'> or <class 'astropy.units.core.CompositeUnit'>)

    Raises
    ------
    ValueError
        Because `unit` did not parse as unit.

    Notes
    ----
    For the complete list of HEK parameters: https://www.lmsal.com/hek/VOEvent_Spec.html

    """
    cm2 = u.def_unit("cm2", u.cm**3)
    m2 = u.def_unit("m2", u.m**2)
    m3 = u.def_unit("m3", u.m**3)

    aliases = {
        "steradian": u.sr,
        "arcseconds": u.arcsec,
        "degrees": u.deg,
        "sec": u.s,
        "emx": u.Mx,
        "amperes": u.A,
        "ergs": u.erg,
        "cubic centimeter": u.ml,
        "square centimeter": cm2,
        "cubic meter": m3,
        "square meter": m2,
    }

    with u.add_enabled_units([cm2, m2, m3]), u.set_enabled_aliases(aliases):
        # If they are units of coordinates, it will have more than one unit,
        # otherwise it will be just one unit.
        # NOTE: There is an assumption that coord1_unit, coord2_unit and coord3_unit will be the same.
        units = re.split(r'[, ]', unit)
        return u.Unit(units[0].lower())

Another thing that has been done was adding a documentation string for parse_chaincode function.

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

    Parameters
    ----------
    value: PolygonSkyRegion
        A polygon defined using vertices in sky coordinates.
    attribute: dict
        An object from coord_properties.json
    unit: str
        The unit of the coordinates

    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.

    """
    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)

    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, representation_type="cylindrical", frame="heliocentric")
    else:
        vertices = SkyCoord(coord1_list, coord2_list, frame=attribute["frame"])

    return PolygonSkyRegion(vertices = vertices)

Test

Ahmed Hossam — Sun, 12 May 2024 20:19:12 +0000

This is a test post to see if RSS works.