Plotting maps from different times on same Carrington map

Hi, I have a project where I want to plot the contours of coronal holes from AIA data over HMI data. The catch is that I’m looking at farside HMI data and I want to compare AIA data that is collected half a Carrington rotation earlier (when the farside was, well, nearside). I have rotated both datasets into the Carrington frame using make_heliographic_header and reproject_to. When I plot them side-by-side, they look like they’re on the same coordinate system. I can overlay the AIA data onto the farside HMI data axis if I set the autoalign keyword to False.

However, I really want the contours from the AIA data, and whenever I try to call draw_contours for the AIA data onto the axis with farside HMI data, they are shifted by ~180 degrees:

I suspect this has to do with the different observation times. Is there any way to overcome this? Thanks in advance! (Sorry for the low information, I wanted to add more pictures but I can’t since I’m a new user.)

I’ve tracked this problem down to a coordinate conversion problem. It appears that if you plot a calculated HGC coordinate onto a map in HGC coordinates with data from a different time, it will recalculate the former HGC coordinate by taking the Stonyhurst reference and rotating it by the new HGC map.

Here is a demonstration of the problem in code:

import astropy.units as u
from astropy.coordinates import SkyCoord
from sunpy.coordinates import frames


# Make Stonyhurst frame at a time:
point_time = '2023-12-07 13:36'
stonyhurst_frame = frames.HeliographicStonyhurst(
    obstime=point_time)
# Create the point represented the launch lon/lat (HEEQ = Stonyhurst)
point_in_stonyhurst = SkyCoord(
    lon=(170 * u.deg), lat=(-1 * u.deg),
    frame=stonyhurst_frame)
print("Stonyhurst coordinate @ Ti:",
      point_in_stonyhurst)

# Get the true Carrington coordinate for this time:
carrington_frame = frames.HeliographicCarrington(
    obstime=point_time, observer='earth')
carrington_point = point_in_stonyhurst.transform_to(carrington_frame)
print("HGC coordinate @ Ti:",
      carrington_point)

# This will ROTATE if you transform this into an HGC frame that
# has a different time!
different_time = '2023-11-23 12:00'
carrington_frame_different_time = frames.HeliographicCarrington(
    obstime=different_time, observer='earth')
carrington_point_different_time =\
    carrington_point.transform_to(carrington_frame_different_time)
print("HGC coordinate @ Ti -> HGC coordinate @ Tf:\n",
      carrington_point_different_time)

If you run this, the coordinate that is rotated in HGC changes if it is rotated in HGC again. This doesn’t make sense to me – a coordinate already in HGC should stay the same in any HGC frame, regardless of observation info.

There is also no way around this rotation occurring when you plot data, e.g. with draw_contours, because it always runs a rotation if the WCS is different, which is always true for Carrington frames at different times because the observation time is different.

Is there a way around this?

I’ll address this specific point first. A coordinate is by default treated as if it indicates a point in inertial space. For example, if a coordinate points to an active region at one time, at a later time it will not continue to point to that active region because the Sun will have rotated that active region away. That’s why a coordinate in HGC changes over time, again, by default.

For your overall issue, you basically have two options:

  1. You can wrap all your plotting code in the propagate_with_solar_surface() context manager, and then (differential) solar rotation will be applied to all coordinate transformations. That is, coordinates on the solar surface will track a solar feature. I’ll note that only if you choose rotation_model='rigid' – which is not the default rotation model – will HGC coordinates not evolve at all over time; otherwise solar differential rotation will cause coordinates to drift.
  2. You can force the WCS information for your two maps to be the same after you have performed the reprojections. The most expedient way to do this is to create a new map by combining the data of one map with the metadata of the other map, e.g.:
aia_map_new = sunpy.map.Map((aia_map.data, hmi_map.meta))

Feel free to ask follow-up questions!

2 Likes