# Derotate HMI map at carrington rate with Postel projection

**URL:** <https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840>\
**Category:** SunPy\
**Tags:** sunpy\
**Created:** [December 6, 2023, 5:32pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840 "2023-12-06T17:32:21Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Wall.E](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/wall.e/32/146_2.png) [@Wall.E](https://community.openastronomy.org/u/Wall.E)\
**Post date:** [December 6, 2023, 5:32pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/1 "2023-12-06T17:32:21Z")

</div>

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?

```python
# 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)

```

---

<div class="post-metadata">

**Author:** ![ayshih](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/ayshih/32/84_2.png) [@ayshih](https://community.openastronomy.org/u/ayshih)\
**Post date:** [December 6, 2023, 5:49pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/2 "2023-12-06T17:49:52Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Wall.E](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/wall.e/32/146_2.png) [@Wall.E](https://community.openastronomy.org/u/Wall.E)\
**Post date:** [December 6, 2023, 5:59pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/3 "2023-12-06T17:59:16Z")

</div>

> [@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?

---

<div class="post-metadata">

**Author:** ![ayshih](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/ayshih/32/84_2.png) [@ayshih](https://community.openastronomy.org/u/ayshih)\
**Post date:** [December 8, 2023, 3:54pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/4 "2023-12-08T15:54:21Z")

</div>

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:

```python
# 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:

```python
# 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.)

---

<div class="post-metadata">

**Author:** ![Wall.E](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/wall.e/32/146_2.png) [@Wall.E](https://community.openastronomy.org/u/Wall.E)\
**Post date:** [December 12, 2023, 6:16am UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/5 "2023-12-12T06:16:22Z")

</div>

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:

```python
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()

```

 ![Screenshot 2023-12-11 at 10.08.58 PM](https://europe1.discourse-cdn.com/flex005/uploads/openastronomy/original/1X/93689f234807b9f80c314ac984821f6ab2cefe69.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![ayshih](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/ayshih/32/84_2.png) [@ayshih](https://community.openastronomy.org/u/ayshih)\
**Post date:** [December 12, 2023, 4:14pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/6 "2023-12-12T16:14:24Z")

</div>

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](https://docs.sunpy.org/en/stable/topic_guide/rsun.html) 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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![Wall.E](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/wall.e/32/146_2.png) [@Wall.E](https://community.openastronomy.org/u/Wall.E)\
**Post date:** [December 12, 2023, 7:03pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/7 "2023-12-12T19:03:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ayshih](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/ayshih/32/84_2.png) [@ayshih](https://community.openastronomy.org/u/ayshih)\
**Post date:** [December 12, 2023, 8:03pm UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/8 "2023-12-12T20:03:14Z")

</div>

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](https://github.com/sunpy/sunpy/issues/7345)) right after my earlier response.

---

<div class="post-metadata">

**Author:** ![Wall.E](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/wall.e/32/146_2.png) [@Wall.E](https://community.openastronomy.org/u/Wall.E)\
**Post date:** [December 13, 2023, 1:22am UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/9 "2023-12-13T01:22:25Z")

</div>

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

```python
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?

```auto
# 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()

```

 ![Screenshot 2023-12-12 at 5.21.06 PM](https://europe1.discourse-cdn.com/flex005/uploads/openastronomy/original/1X/9593ee4c808213d94d1ac5fb7a877b88de76af80.jpeg)

---

<div class="post-metadata">

**Author:** ![ayshih](https://dub1.discourse-cdn.com/flex005/user_avatar/community.openastronomy.org/ayshih/32/84_2.png) [@ayshih](https://community.openastronomy.org/u/ayshih)\
**Post date:** [December 13, 2023, 2:22am UTC](https://community.openastronomy.org/t/derotate-hmi-map-at-carrington-rate-with-postel-projection/840/10 "2023-12-13T02:22:03Z")

</div>

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