Derotate HMI map at carrington rate with Postel projection

[Wall.E]

I would like to derotate a series of HMI images at the Carrington rate, and have them reprojected with the Postel projection at a custom origin. I do not want to derotate differentially; so as to keep the relative effects of differential rotation in my reprojected series of images.
The documentation is explicit about differentially derotating, but I’m not seeing anything for a more “rigid” derotation at a constant Carrington rate.

Should I simply build a header with an origin given in Carrington coordinates that remains constant for all my maps when using reproject_to()?

For example, looping what’s below over all my images?

# Define origin in hpc first 
origin_hpc = SkyCoord(735*u.arcsec, -340*u.arcsec, frame=hmi_map.coordinate_frame)
# Convert it to Carrington coordinates
origin = origin_hpc.heliographic_carrington

out_shape = (750, 750)
out_header = sunpy.map.make_fitswcs_header(
    out_shape,
    origin,
    scale=[0.03, 0.03]*u.deg/u.pix,
    projection_code="ARC"
)

out_map = hmi_map.reproject_to(out_header)

[ayshih]

You didn’t specify the approach you are using to derotate the images, but both the propagate_with_solar_surface() context manager and the RotatedSunFrame metaframe class support the keyword argument rotation_model='rigid' to specify rigid (de)rotation. Does that accomplish what you want?

[Wall.E]

@ayshih

You didn’t specify the approach you are using to derotate the images

That’s what my approach above is asking: if all the reprojected maps are centered at a constant Carrington coordinate, wouldn’t they by definition derotated at Carrington rate?

With your approach using propagate_with_solar_surface(), I guess the Postel projection will be a second step, but then again, I would have to give it an origin that remains constant in Carrington coordinates, so there will be yet again a reprojection, I feel like the job would be done twice, wouldn’t it?

[ayshih]

Yes, broadly speaking there are two ways to approach this.

Define a Carrington WCS header for each map

This is what you would do to loop over the code snippet you provided:

# Choose a point in Carrington heliographic coordinates using the first map
origin_hpc = SkyCoord(735*u.arcsec, -340*u.arcsec, frame=list_of_maps[0].coordinate_frame)
origin_hgc = origin_hpc.heliographic_carrington

out_shape = (750, 750)

for this_map in list_of_maps:
    # Define the map-specific origin using the above Carrington heliographic coordinate
    this_origin = SkyCoord(origin_hgc.lon, origin_hgc.lat, frame='heliographic_carrington',
                           obstime=this_map.date, observer=this_map.observer_coordinate)
    out_header = sunpy.map.make_fitswcs_header(
        out_shape,
        this_origin,
        scale=[0.03, 0.03]*u.deg/u.pix,
        projection_code="ARC"
    )

    out_map = this_map.reproject_to(out_header)
    
    plt.figure()
    out_map.plot()

You’d first calculate the Carrington longitude/latitude of interest, and then for each map you’d create a map-specific WCS header. Note that while each reprojected map will look like how you want, they are formally each in a different coordinate frame. You may then find yourself wanting to access the data array directly (.data) to avoid some of the “smarts” in sunpy.

Define a single Carrington WCS header

Alternatively, you can leverage propagate_with_solar_surface(rotation_model='rigid') to get reprojected maps in the same coordinate frame by using a single Carrington WCS header:

# Choose a point in Carrington heliographic coordinates using the first map
origin_hpc = SkyCoord(735*u.arcsec, -340*u.arcsec, frame=list_of_maps[0].coordinate_frame)
origin_hgc = origin_hpc.heliographic_carrington

out_shape = (750, 750)
out_header = sunpy.map.make_fitswcs_header(
    out_shape,
    origin_hgc,
    scale=[0.03, 0.03]*u.deg/u.pix,
    projection_code="ARC"
)

for this_map in list_of_maps:
    with propagate_with_solar_surface(rotation_model='rigid'):
        out_map = this_map.reproject_to(out_header)
    
    plt.figure()
    out_map.plot()

Now the output maps are all in the same coordinate frame: Carrington heliographic coordinates for the observing time and observer location of the first map.

(Essentially, the use of propagate_with_solar_surface(rotation_model='rigid') can be thought of as telling sunpy that the desired transformation of Carrington coordinates at one time to Carrington coordinates at another time should be the simplest one where the longitude/latitude values don’t change.)

[Wall.E]

That’s what I’ve been trying (your first method), and there seems to be a problem with the reprojection… It removes data above the origin point, where there should be data (see image down below). I checked with the reprojection pipeline from JSOC with the same origin and shape parameters coordinates, and I don’t have this problem.

Here’s a reproducible example:

import matplotlib
import matplotlib.pyplot as plt
import sunpy.map
import astropy.units as u
from astropy.coordinates import SkyCoord
from sunpy.coordinates import frames
import sunpy.data.sample

hmi_map = sunpy.map.Map(sunpy.data.sample.HMI_LOS_IMAGE)
hmi_rotated = hmi_map.rotate(order=3)
# Define center of reprojected map
origin = SkyCoord(0*u.deg, 50*u.deg, frame=frames.HeliographicStonyhurst, obstime=hmi_map.date, 
                  observer=hmi_rotated.observer_coordinate)

hmi_map = sunpy.map.Map(sunpy.data.sample.HMI_LOS_IMAGE)
hmi_rotated = hmi_map.rotate(order=3)
# Define center of reprojected map
origin = SkyCoord(0*u.deg, 50*u.deg, frame=frames.HeliographicStonyhurst, obstime=hmi_map.date, 
                  observer=hmi_rotated.observer_coordinate)
# Define header for the reprojection
out_shape = (1024, 1024)
out_header = sunpy.map.make_fitswcs_header(
    out_shape,
    origin,
    scale=[0.0301, 0.0301]*u.deg/u.pix,
    projection_code="ARC" #  azimuthal equidistant projection, aka Postel
    )
# Reproject
out_map = hmi_map.reproject_to(out_header)

fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 2, 1, projection=hmi_rotated)
hmi_rotated.plot(axes=ax)
hmi_rotated.draw_grid(axes=ax, color='blue')
hmi_rotated.draw_limb(axes=ax, color='blue')
ax.plot_coord(origin, 'o', color='red', fillstyle='none', markersize=20)

ax = fig.add_subplot(1, 2, 2, projection=out_map)
out_map.plot(axes=ax)
out_map.draw_grid(axes=ax, color='blue')
out_map.draw_limb(axes=ax, color='blue')
ax.plot_coord(origin, 'o', color='red', fillstyle='none', markersize=20)
ax.set_title('Postel projection centered at ROI')
plt.show()

As I’m using a full-disk image (HMI magnetogram), we have data all the way near the pole, whereas on the right, the reprojected map has removed the data much below that, with NaN values, around 45 degrees, as if it was putting the north limb down there. So, there must be something wrong with reproject(). @svank That’s what I mentionned to you at AGU.

[ayshih]

You’re yet another person that’s been bitten by the fact that SDO uses a different radius for the Sun than the default value in sunpy (696 Mm versus 695.7 Mm). Normally this mismatch manifests as “lost” data near the limb, but the pixels in your reprojected map are so small that you’re losing data even pretty far from the limb. See The role of “rsun” in sunpy — sunpy 5.1.0 documentation for relevant explanation.

In your case, the solution is to explicitly set the rsun frame attribute of origin (i.e., the reference coordinate of your output map) to match the rsun of the original coordinate frame:

origin = SkyCoord(0*u.deg, 50*u.deg, frame=frames.HeliographicStonyhurst, obstime=hmi_map.date, 
                  observer=hmi_rotated.observer_coordinate, rsun=hmi_rotated.coordinate_frame.rsun)

All of the data will survive the reprojection.

[Wall.E]

Many thanks, that worked.
Wouldn’t it better to have reproject_to() issue a warning at least for this problem? Reproject has the map object, and the origin. If I don’t give rsun in when defining the origin with SkyCoord, the Sunpy default is assumed. But reproject_to() is given a map, which as the rsun attribute from HMI, so reproject could check and warn of any discrepancy between rsun in origin and the one in the map object.

[ayshih]

Heh, in fact I created a GitHub issue for precisely that (Show a warning when calling `reproject_to()` and there is a (minor) mismatch in `rsun` · Issue #7345 · sunpy/sunpy · GitHub) right after my earlier response.

[Wall.E]

After I have done some processing on the data of the reprojected Postel map, for example:

import matplotlib
import matplotlib.pyplot as plt
import sunpy.map
import astropy.units as u
from astropy.coordinates import SkyCoord
from sunpy.coordinates import frames
import sunpy.data.sample

hmi_map = sunpy.map.Map(sunpy.data.sample.HMI_LOS_IMAGE)
hmi_rotated = hmi_map.rotate(order=3)
# Define center of reprojected map
origin = SkyCoord(0*u.deg, 50*u.deg, frame=frames.HeliographicStonyhurst, obstime=hmi_map.date, 
                  observer=hmi_rotated.observer_coordinate, rsun=hmi_map.coordinate_frame.rsun)
# Define header for the reprojection
out_shape = (1024, 1024)
out_header = sunpy.map.make_fitswcs_header(
    out_shape,
    origin,
    scale=[0.0301, 0.0301]*u.deg/u.pix,
    projection_code="ARC" #  azimuthal equidistant projection, aka Postel
    )
# Reproject
out_map = hmi_map.reproject_to(out_header)
# Do something with the reprojected patch
out_map2 = out_map**3

if I want to reproject again that processed data onto the original Helioprojective grid of the original full-disk magnetogram, should I just use reproject_to() using the wcs attribute of the original full-disk map from which the postel patch was made, like this?

# Reproject on whole solar disk
out_map_disk = out_map2.reproject_to(hmi_rotated.wcs)

ms=6

# Plot
fig = plt.figure(figsize=(18, 7))
ax1 = fig.add_subplot(1, 3, 1, projection=hmi_rotated)
hmi_rotated.plot(axes=ax1)
hmi_rotated.draw_grid(axes=ax1, color='blue')
hmi_rotated.draw_limb(axes=ax1, color='blue')
ax1.plot_coord(origin, 'o', color='red', fillstyle='none', markersize=ms)

ax2 = fig.add_subplot(1, 3, 2, projection=out_map)
out_map.plot(axes=ax2)
out_map.draw_grid(axes=ax2, color='blue')
out_map.draw_limb(axes=ax2, color='blue')
ax2.plot_coord(origin, 'o', color='red', fillstyle='none', markersize=ms)


ax3 = fig.add_subplot(1, 3, 3, projection=out_map_disk)
out_map_disk.plot(axes=ax3)
out_map_disk.draw_grid(axes=ax3, color='blue')
out_map_disk.draw_limb(axes=ax3, color='blue')
ax3.plot_coord(origin, 'o', color='red', fillstyle='none', markersize=ms)

plt.show()

[ayshih]

Yup. Keep in mind for your processing in Postel space that the (default) reprojection back and forth is merely an interpolation.