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_coord1andevent_coord2andevent_coord3into oneSkyCoordobject 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.
Top comments (0)