Difference between propagate_with_solar_surface and transform_with_sun_center

[GiuliaMRonca]

Hello,

I want to differentially rotate an image, using the procedure explained in Differentially rotating a map — SunPy 4.1.0 documentation

I used to do that one year ago with this procedure:

# Reproject the map from the input frame to the output frame.
with transform_with_sun_center():
    arr, _ = reproject_interp(mapname, out_wcs, out_shape)
# Create the output map and preserve the original map’s plot settings.
out_warp = sunpy.map.Map(arr, out_wcs)

This is what was suggested at the time. Now, I see that the suggested function is:

with propagate_with_solar_surface():
out_warp = aiamap.reproject_to([out_wcs])

I have read about the difference in the dedicated pages, but it is still not clear to me what exactly changes in one case or the other. Does this upgrade mean that what I did previously with the transform_with_sun_center function is wrong?

[ayshih]

Thanks for asking!

First, I’ll clarify that the transform_with_sun_center() context manager is solely about ignoring the slow translational motion of the Sun through inertial space. (The propagate_with_solar_surface() context manager ignores this translational motion as well.)
transform_with_sun_center() by itself does not incorporate any notion of differential rotation.

The “old” approach for differentially rotating a map actually requires an important modification to the WCS instance prior to your quoted code snippet:

out_wcs = WCS(...)

# Assuming appropriate definitions of out_frame and in_time
from sunpy.coordinates import RotatedSunFrame
rot_frame = RotatedSunFrame(base=out_frame, rotated_time=in_time)
out_wcs.coordinate_frame = rot_frame

with transform_with_sun_center():
    arr, _ = reproject_interp(aiamap, out_wcs, out_shape)
out_warp = sunpy.map.Map(arr, out_wcs)

and should be contrasted against the “new” approach:

out_wcs = WCS(...)

# No modification to out_wcs is needed

with propagate_with_solar_surface():
    arr, _ = reproject_interp(aiamap, out_wcs, out_shape)
out_warp = sunpy.map.Map(arr, out_wcs)

which is further even simpler using the now-available reproject_to() method:

out_wcs = WCS(...)

with propagate_with_solar_surface():
    out_warp = aiamap.reproject_to(out_wcs)

In other words:

• With the propagate_with_solar_surface() context manager, it is no longer necessary to modify the WCS instance to include a RotatedSunFrame instance.
• With the reproject_to() method, it is no longer necessary to manually call reproject_interp() and handle the returned data array.

None of the reprojected output has changed, and your old code is likely not wrong (but feel free to share it so that it can be checked!). Switching to the current approach is simply about simplifying your code.

[GiuliaMRonca]

Hi @ayshih ,
thanks for your reply. It was really helpful.

I also tried to subtract one image from the other to qualitatively check the difference in the approaches.
I am going to share my code. If you could tell me your opinion, I would love to read it.

I have a set of sunpy maps (here after called
raw=sunpy.map.Map(event+‘/raw_img/*.fits’, sequence=True)
). They are images taken on March 07 2012.
raw[0] corresponds to 2012-03-06 at 23:32:00.620
while raw[i=350] corresponds to 2012-03-09 at 11:32:00.62

base_raw=raw[0]
base_updated= update_pointing(base_raw)
base_registered = register(base_updated)
basemap= normalize_exposure(base_registered)
new_dimensions = [2048, 2048] * u.pixel # Resampling
basemap = basemap.resample(new_dimensions)

i=350

if raw[i].meta[‘exptime’]>2.9 and raw[i].meta[‘exptime’]<3:

smap_raw =raw[i]
smap_updated_pointing_old= update_pointing(smap_raw)
smap_registered = register(smap_updated_pointing_old)
smap= normalize_exposure(smap_registered)
new_dimensions = [2048, 2048] * u.pixel
map_old = smap.resample(new_dimensions)

Then I differentially rotate smap_raw first with both the methods:

base_time= sunpy.time.parse_time(basemap.meta[‘date-obs’])
t= sunpy.time.parse_time(map_old.meta[‘date-obs’])
time_diff= -(t-base_time).to(u.min)
out_old=diff_rot(map_old, time_diff, basemap)
out_new=diff_rot_new(map_old, time_diff, basemap)

The two functions are below:

def diff_rot(mapname, time_diff, basemap):

#NOTE: to go back in time, time_diff MUST be negative
in_time = mapname.date
out_time = in_time + time_diff
out_frame = Helioprojective(observer='earth', obstime= out_time, rsun= mapname.meta['rsun_ref']*u.m)
rot_frame = RotatedSunFrame(base=out_frame, rotated_time=in_time)

out_shape = mapname.data.shape
# WCS converts cdelt and cunit from arcsec to deg
out_wcs = basemap.wcs
out_wcs.coordinate_frame = rot_frame

# Reproject the map from the input frame to the output frame.
with transform_with_sun_center():
    arr, _ = reproject_interp(mapname, out_wcs, out_shape)
# Create the output map and preserve the original map’s plot settings.
out_warp = sunpy.map.Map(arr, out_wcs)
out_warp.plot_settings = mapname.plot_settings


mapname.meta["date-obs"]=mapname.meta["date-obs"]
if mapname.name[0:3]=='AIA':
    out_warp.meta["waveunit"] = mapname.meta["waveunit"]
    out_warp.meta["exptime"] = mapname.meta['exptime']
    out_warp.meta["wavelnth"] = mapname.meta['wavelnth']
    out_warp.meta["instrume"]= mapname.meta['instrume']
    out_warp.meta["telescop"]= mapname.meta['telescop']
    out_warp.meta["rsun_obs"]= mapname.meta['rsun_obs']

out_warp.meta["cdelt1"]=mapname.meta['cdelt1']
out_warp.meta["cdelt2"]=mapname.meta['cdelt2']
out_warp.meta["cunit1"]=mapname.meta['cunit1']
out_warp.meta["cunit2"]=mapname.meta['cunit2']

return out_warp

And here is the “new” version of the differential rotation function:

from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

import sunpy.map
from sunpy.coordinates import Helioprojective, propagate_with_solar_surface

def diff_rot_new(mapname, time_diff, basemap):

#NOTE: to go back in time, time_diff MUST be negative
in_time = mapname.date
out_time = in_time + time_diff

out_frame = Helioprojective(observer='earth', obstime= out_time, rsun= mapname.meta['rsun_ref']*u.m)


# WCS converts cdelt and cunit from arcsec to deg
out_wcs = basemap.wcs

# Reproject the map from the input frame to the output frame.
with propagate_with_solar_surface():
    out_warp = mapname.reproject_to(out_wcs)

out_warp.meta["date-obs"]=mapname.meta["date-obs"]
if mapname.name[0:3]=='AIA':
    out_warp.meta["waveunit"] = mapname.meta["waveunit"]
    out_warp.meta["exptime"] = mapname.meta['exptime']
    out_warp.meta["wavelnth"] = mapname.meta['wavelnth']
    out_warp.meta["instrume"]= mapname.meta['instrume']
    out_warp.meta["telescop"]= mapname.meta['telescop']
    out_warp.meta["rsun_obs"]= mapname.meta['rsun_obs']
    out_warp.meta["pc1_1"]= mapname.meta['pc1_1']
    out_warp.meta["pc1_2"]= mapname.meta['pc1_2']
    out_warp.meta["pc2_1"]= mapname.meta['pc2_1']
    out_warp.meta["pc2_2"]= mapname.meta['pc2_2']
    
out_warp.meta["cdelt1"]=mapname.meta['cdelt1']
out_warp.meta["cdelt2"]=mapname.meta['cdelt2']
out_warp.meta["cunit1"]=mapname.meta['cunit1']
out_warp.meta["cunit2"]=mapname.meta['cunit2']

return out_warp

As I told you, I then tried to subtract one image from the other:

diff = out_new.data - out_old.data

You can see the results.
There is some misalignement in the pixel positions so the difference’s value is quite high sometimes, especially close to the active regions, where bright and dark pixels are very close to each other.
Here I limited the plot to -7 and +7 values (maximum value is +155, minimum value is -193)

[GiuliaMRonca]

In percentage, the difference is interesting, too.

diff = out_new.data - out_old.data
perc= ( diff / out_new.data )*100
meta = out_new.meta
perc_map = sunpy.map.Map(perc, meta)

I think this is mostly due to a misalignement of the pixels, rather than an error, but coul this lead to any problem in the future? Especially when different people process images with different methods

[ayshih]

The reason that there is a discrepancy between the outputs from your two functions is because there is a difference in the target coordinate frame. In short, diff_rot() reprojects to an observer at Earth center, while diff_rot_new() reprojects to an observer at SDO’s location. Accordingly, the reprojection is slightly different.

diff_rot()

You have:

out_frame = Helioprojective(observer='earth', obstime= out_time, rsun= mapname.meta['rsun_ref']*u.m)

Since out_time = in_time + time_diff = mapname.date + (basemap.date - mapname.date) = basemap.date, we can rewrite this as:

out_frame = Helioprojective(observer='earth', obstime=basemap.date, rsun=basemap.coordinate_frame.rsun)

However, note that the observer is specified as Earth and not SDO. The analogous call if you want the observer to be SDO would be:

out_frame = Helioprojective(observer=basemap.observer_coordinate, obstime=basemap.date, rsun=basemap.coordinate_frame.rsun)

or far more simply:

out_frame = basemap.coordinate_frame

Because you explicitly attach a coordinate frame to out_wcs, the reprojection will use that coordinate frame – with the observer at Earth – instead of accessing the observer information contained in out_wcs.

diff_rot_new()

You have out_frame defined the same way, with the observer specified as Earth instead of as SDO. However, defining out_frame is misleading because it is never used. Instead, because you reproject directly to out_wcs, which contains the observer location as SDO, the target coordinate frame has the observer as SDO.

[GiuliaMRonca]

I see.
It is perfect now and my doubts are cleared.

Thank you very much!