Coordinate frames that include solar differential rotation #3537

ayshih · 2019-11-22T19:08:21Z

This PR integrates solar (differential) rotation into the coordinates framework. The approach here is to enable defining a frame that is a modified version of a base coordinate frame by a duration of solar (differential) rotation. The class factory RotatedSunFrame is supplied the "base" coordinate frame, and creates a base-specific rotated-Sun frame class.

Closes #3357

Still needs:

Add tutorial (see this recent build)
Finish writing tests
Support array obstime/duration/rotated_time

Walkthrough

>>> import astropy.units as u
>>> from sunpy.coordinates import RotatedSunFrame
>>> import sunpy.coordinates.frames as f

Propagating a coordinate

Let's define a coordinate array of points in a rotated-Sun frame that in the base HGS frame were on the meridian as seen from Earth. Here, the rotated-Sun frame applies +1 day of solar rotation.

>>> hgs = f.HeliographicStonyhurst(obstime='2001-01-01')
>>> rotated_arc = RotatedSunFrame([0]*11*u.deg, range(-75, 90, 15)*u.deg, base=hgs, duration=1*u.day)
>>> rotated_arc
<RotatedSunHeliographicStonyhurst Coordinate (base=<HeliographicStonyhurst Frame (obstime=2001-01-01T00:00:00.000)>, duration=1.0 d, rotation_model=howard): (lon, lat, radius) in (deg, deg, km)
    [(0., -75., 695700.), (0., -60., 695700.), (0., -45., 695700.),
     (0., -30., 695700.), (0., -15., 695700.), (0.,   0., 695700.),
     (0.,  15., 695700.), (0.,  30., 695700.), (0.,  45., 695700.),
     (0.,  60., 695700.), (0.,  75., 695700.)]>

The coordinate components are relative to the distorted axes of the rotated-Sun frame, and thus are already differentially rotated in real space. Convert the coordinate array back to the HGS frame to see what I mean.

>>> rotated_arc.transform_to(hgs)
<HeliographicStonyhurst Coordinate (obstime=2001-01-01T00:00:00.000): (lon, lat, radius) in (deg, deg, km)
    [(10.7550473 , -75., 695700.), (11.70697161, -60., 695700.),
     (12.80904447, -45., 695700.), (13.68216339, -30., 695700.),
     (14.17617983, -15., 695700.), (14.32632838,   0., 695700.),
     (14.17617983,  15., 695700.), (13.68216339,  30., 695700.),
     (12.80904447,  45., 695700.), (11.70697161,  60., 695700.),
     (10.7550473 ,  75., 695700.)]>

Another example

There are several ways to instantiate a RotatedSunFrame frame. Here, the input supplied to base has the relevant coordinate data, so that data is automatically "moved" to the correct place in the resulting frame. Also, instead of supplying duration directly, the rotated_time keyword is used to indirectly define the duration of rotation.

>>> center = f.Helioprojective(0*u.arcsec, 0*u.arcsec, obstime='2001-01-01', observer='earth')
>>> rotated_center = RotatedSunFrame(base=center, rotated_time='2001-01-03')
>>> rotated_center
<RotatedSunHelioprojective Coordinate (base=<Helioprojective Frame (obstime=2001-01-01T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>)>, duration=2.0 d, rotation_model=howard): (Tx, Ty) in arcsec
    (0., 0.)>

No surprises here: The HPC center of the solar disk moves mostly in the X direction after +2 days of solar rotation.

>>> rotated_center.transform_to(center)
<Helioprojective Coordinate (obstime=2001-01-01T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>): (Tx, Ty, distance) in (arcsec, arcsec, AU)
    (468.86076705, -6.35238293, 0.97923501)>

Because RotatedSunFrame frames are in the transformation graph, you can do more creative transformations. Here, we calculate what HPC coordinate when differentially rotated by -1 days corresponds to the disk center when differentially rotated by +2 days.

>>> rotated_center.transform_to(RotatedSunFrame(base=center, duration=-1*u.day))
<RotatedSunHelioprojective Coordinate (base=<Helioprojective Frame (obstime=2001-01-01T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>)>, duration=-1.0 d, rotation_model=howard): (Tx, Ty, distance) in (arcsec, arcsec, AU)
    (666.1751848, -13.91351551, 0.97991384)>

Grid overlay on a Map

>>> from sunpy.map import Map
>>> from sunpy.data.sample import AIA_171_IMAGE
>>> m = Map(AIA_171_IMAGE)

>>> import matplotlib.pyplot as plt
>>> ax = plt.subplot(projection=m)
>>> m.plot()

>>> overlay = ax.get_coords_overlay('heliographic_stonyhurst')
>>> overlay[0].set_ticks(spacing=15. * u.deg)
>>> overlay[1].set_ticks(spacing=90. * u.deg)
>>> overlay.grid(ls='-', color='white')

>>> overlay = ax.get_coords_overlay(RotatedSunFrame(base=f.HeliographicStonyhurst(obstime=m.date), duration=27*u.day))
>>> overlay[0].set_ticks(spacing=15. * u.deg)
>>> overlay[1].set_ticks(spacing=90. * u.deg)
>>> overlay.grid(ls='-', color='blue')

pep8speaks · 2019-11-22T19:08:26Z

Hello @ayshih! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

In the file sunpy/coordinates/tests/strategies.py:

Line 25:9: E128 continuation line under-indented for visual indent

In the file sunpy/coordinates/tests/test_metaframes.py:

Line 102:20: E721 do not compare types, use 'isinstance()'
Line 163:20: E721 do not compare types, use 'isinstance()'

Comment last updated at 2020-02-05 15:01:21 UTC

ghost · 2019-11-22T19:08:26Z

Thanks for the pull request @ayshih! Everything looks great!

ayshih · 2019-11-22T21:19:02Z

Multiple rotation durations

You can work with array obstime/duration/rotated_time.

>>> hpc = Helioprojective(obstime='2001-01-01', observer='earth')
>>> sc = SkyCoord([0]*9*u.arcsec, [0]*9*u.arcsec, frame=hpc)
>>> series = RotatedSunFrame(base=sc, duration=range(0, 27, 3)*u.day)
>>> series
<RotatedSunHelioprojective Coordinate (base=<Helioprojective Frame (obstime=2001-01-01T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>)>, duration=[ 0.  3.  6.  9. 12. 15. 18. 21. 24.] d, rotation_model=howard): (Tx, Ty) in arcsec
    [(0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
     (0., 0.), (0., 0.)]>
>>> series.transform(hpc)
<Helioprojective Coordinate (obstime=2001-01-01T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'earth'>): (Tx, Ty, distance) in (arcsec, arcsec, AU)
    [(   0.        ,    0.        , 0.97866506),
     ( 666.1751848 ,  -13.91351551, 0.97991384),
     ( 971.99406549,  -48.04083807, 0.98298359),
     ( 756.05760587,  -83.90802534, 0.98622017),
     ( 137.55563461, -102.40288732, 0.98789324),
     (-553.88198072,  -93.79075537, 0.98711382),
     (-951.06271076,  -62.59442403, 0.9842956 ),
     (-839.39922086,  -25.34560595, 0.98094109),
     (-274.84399172,   -2.09057135, 0.9788526 )]>

ayshih · 2019-11-23T00:30:59Z

Transforming a coordinate to a frame with a different observation time

In the coordinates framework, a coordinate (e.g., SkyCoord) points to a specific location in inertial space, while the obstime/observer define the origin and axis orientations of the coordinate frame. If you transform a SkyCoord to a coordinate frame with a different obstime, you obtain the same inertial location in the new coordinate frame, because there there isn't any way to specify how that coordinate evolves in time.

For example, here we try to transform disk center as seen by one observer to an observer 90 degrees separated in longitude and 6 days later in time:

>>> observer1 = SkyCoord(HeliocentricInertial(0*u.deg, 0*u.deg, 1*u.AU, obstime='2001-01-01'))
>>> observer2 = SkyCoord(HeliocentricInertial(90*u.deg, 0*u.deg, 1*u.AU, obstime='2001-01-07'))

>>> frame1 = Helioprojective(observer=observer1, obstime=observer1.obstime)
>>> frame2 = Helioprojective(observer=observer2, obstime=observer2.obstime)

>>> sc = SkyCoord(0*u.arcsec, 0*u.arcsec, frame=frame1)
>>> sc.transform_to(frame2)
<SkyCoord (Helioprojective: obstime=2001-01-07T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2001-01-07T00:00:00.000): (lon, lat, radius) in (deg, deg, AU)
    (59.22396701, 0., 1.)>): (Tx, Ty, distance) in (arcsec, arcsec, AU)
    (-956.93749939, 1.16033026, 1.00006359)>

The output coordinate is on the east limb, which is consistent with not evolving the location with solar rotation.

With this PR, you can now include the evolution of solar rotation via this call:

>>> RotatedSunFrame(base=sc, rotated_time=frame2.obstime).transform_to(frame2)
<Helioprojective Coordinate (obstime=2001-01-07T00:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2001-01-07T00:00:00.000): (lon, lat, radius) in (deg, deg, AU)
    (59.22396701, 0., 1.)>): (Tx, Ty, distance) in (arcsec, arcsec, AU)
    (-78.69026983, 1.16579504, 0.99537569)>

The RotatedSunFrame applies the differential rotation of the 6 days, and then the coordinate is transformed to the new coordinate frame.

sunpy/coordinates/metaframes.py

ayshih · 2019-11-27T18:01:59Z

Differentially rotating a map (versus the old approach)

We can use this frame, plus reproject, to mimic the old approach to produce a differentially rotated map (differential_rotate). This requires reproject 0.6 as well as PR #3540. Importantly, the target frame for the reprojection needs to be the output frame differentially rotated to the time of the input frame.

import matplotlib.pyplot as plt

from reproject import reproject_interp

from astropy.coordinates import SkyCoord
import astropy.units as u
import sunpy.map
from astropy.wcs import WCS
from sunpy.coordinates import RotatedSunFrame, get_body_heliographic_stonyhurst, frames as f
from sunpy.data.sample import AIA_171_IMAGE
from sunpy.physics.differential_rotation import differential_rotate

delta_t = 5*u.day

m = sunpy.map.Map(AIA_171_IMAGE)

# Define the output frame with an associated new observer (at Earth)
output_frame = f.Helioprojective(observer='earth', obstime=m.date + delta_t)

# Define the RotatedSunFrame (its `duration` will turn out to be -1 * `delta_t`)
rot_frame = RotatedSunFrame(base=output_frame, rotated_time=m.date)

out_shape = (1024, 1024)
header = sunpy.map.make_fitswcs_header(out_shape,
                                       m.center.transform_to(rot_frame).as_base(),
                                       scale=u.Quantity(m.scale))

out_wcs = WCS(header)
out_wcs.coordinate_frame = rot_frame

with sunpy.coordinates.transformations.ignore_sun_motion():
    arr, foot = reproject_interp(m, out_wcs, out_shape)

out_warp = sunpy.map.Map(arr, out_wcs)
out_warp.plot_settings = m.plot_settings

old = differential_rotate(m, time=delta_t)

fig = plt.figure(figsize=(15,8))
ax = plt.subplot(1,2,1, projection=out_warp)
out_warp.plot(ax)
ax.set_title("Reproject using RotatedSunFrame")
overlay = ax.get_coords_overlay(
    RotatedSunFrame(base=f.HeliographicStonyhurst(obstime=m.date), rotated_time=output_frame.obstime))
overlay[0].set_ticks(spacing=15. * u.deg)
overlay[1].set_ticks(spacing=90. * u.deg)
with sunpy.coordinates.transformations.ignore_sun_motion():
    overlay.grid(ls='-', color='blue')

ax2 = plt.subplot(1,2,2,projection=old)
old.plot(ax2)
ax2.set_title("Using differential_rotate")

plt.show()

ayshih · 2019-12-18T16:59:14Z

Now rebased off of #3540 so that I can use the context manager in examples. Presumably #3540 will be merged before this PR.

ayshih · 2019-12-18T20:28:24Z

I've added a tutorial (see this recent build).

Cadair · 2019-12-19T12:28:10Z

The tutorial is awesome! Are we confident in the reproject example doing the right thing?

Cadair · 2019-12-19T12:29:03Z

We should also open an issue on astropy about specifying the frame for a WCS. I spoke to both Tom and Erik about it at the astropy meeting and they both agreed it would be something useful to workout upstream.

ayshih · 2019-12-19T15:24:37Z

The tutorial is awesome!

Thanks! I just added another example to show using plot_coord (see recent build).

Are we confident in the reproject example doing the right thing?

Well, I'm confident that the coordinate transformations are correct. Are you confident that reproject is doing the right thing? =P Seriously though, the lingering discrepancies between the "new" method and the "old" method may be difficult to reconcile since we don't actually know that the "old" method was accurate.

ayshih · 2020-01-17T16:08:14Z

Okay, I've committed the test. Who's brave enough to review this? =P

Cadair · 2020-01-17T16:38:49Z

Well it just took me 20 mins to read through your docs so...

Cadair

This is amazing! So cool! I particularly love the wcsaxes example, but I think it is exceptionally well documented and tested. Thanks @ayshih 🚀

Super minimal changes below, but looks good to me. Now, are we going to find a second reviewer for this?!

sunpy/coordinates/metaframes.py

Cadair · 2020-01-17T17:57:17Z

sunpy/coordinates/wcs_utils.py

@@ -15,6 +15,9 @@ def solar_wcs_frame_mapping(wcs):
    type values in the `astropy.wcs.utils.wcs_to_celestial_frame` registry.
    """

+    if hasattr(wcs, "coordinate_frame"):


We really need to upstream this into Astropy...

sunpy/coordinates/tests/test_metaframes.py

Cadair · 2020-01-17T18:06:25Z

Now needs a rebase because #3540 has been merged.

…s coordinate frame attributes

… solar rotation

…otated_time property

dstansby

I haven't thoroughly reviewed the code, but have left some documentation suggestions (and one question).

sunpy/coordinates/frames.py

dstansby · 2020-02-05T15:31:49Z

docs/code_ref/coordinates/rotatedsunframe.rst

@@ -0,0 +1,390 @@
+.. _sunpy-coordinates-rotatedsunframe:
+
+Differential rotation using coordinate frames


I think in this section there should be a clear reference to the available methods (ie. formulae) that can be used for the rotation.

dstansby · 2020-02-05T15:34:28Z

docs/code_ref/coordinates/rotatedsunframe.rst

+
+Creating coordinates
+^^^^^^^^^^^^^^^^^^^^
+


There are lots of class links here that should be intersphinx links and not double-dashes. Happy for another issue to be opened to fix this though.

dstansby · 2020-02-05T15:35:28Z

docs/code_ref/coordinates/rotatedsunframe.rst

+  <HeliographicCarrington Coordinate (obstime=2020-03-04T00:00:00.000): (lon, lat, radius) in (deg, deg, km)
+      (10., 20., 695700.)>
+
+Advanced uses of RotatedSunFrame


I think these examples would be better in individual gallery examples, since this is already quite an individually long bit of documentation.

dstansby

Unless @ayshih wants to put my doc suggestions in here, I'm happy for them to be farmed out to another issue and then another PR.

Cadair · 2020-02-05T18:04:33Z

@dstansby or @ayshih can you open an issue to track the documentation improvements please?

dstansby · 2020-02-05T19:14:13Z

👍 will do

ayshih requested a review from Cadair November 22, 2019 19:08

ayshih added the coordinates Affects the coordinates submodule label Nov 22, 2019

ayshih force-pushed the diffrotdynamic branch 3 times, most recently from dc01ef7 to f7bd66e Compare November 22, 2019 21:18

ayshih force-pushed the diffrotdynamic branch from f7bd66e to a8c3287 Compare November 22, 2019 21:22

Cadair added this to the 2.0 milestone Nov 23, 2019

Cadair reviewed Nov 24, 2019

View reviewed changes

sunpy/coordinates/metaframes.py Outdated Show resolved Hide resolved

sunpy/coordinates/metaframes.py Outdated Show resolved Hide resolved

sunpy/coordinates/metaframes.py Show resolved Hide resolved

sunpy/coordinates/metaframes.py Outdated Show resolved Hide resolved

ayshih force-pushed the diffrotdynamic branch from 57e1fb2 to 21e1c68 Compare November 26, 2019 15:06

wtbarnes mentioned this pull request Dec 17, 2019

MapSequence has a non-functioning keyword derotate #3507

Closed

ayshih force-pushed the diffrotdynamic branch from 21e1c68 to 27754fb Compare December 18, 2019 16:57

ayshih force-pushed the diffrotdynamic branch 3 times, most recently from 81edbcf to e57f274 Compare December 18, 2019 18:58

ayshih force-pushed the diffrotdynamic branch from e57f274 to 16e2a66 Compare December 18, 2019 21:08

ayshih force-pushed the diffrotdynamic branch from 16e2a66 to 784eea8 Compare December 19, 2019 15:07

ayshih force-pushed the diffrotdynamic branch 4 times, most recently from b0257a2 to c0e70d2 Compare December 20, 2019 02:39

ayshih force-pushed the diffrotdynamic branch 2 times, most recently from c35843a to 18fa69f Compare January 17, 2020 16:35

Cadair approved these changes Jan 17, 2020

View reviewed changes

ayshih force-pushed the diffrotdynamic branch from 18fa69f to f103660 Compare January 21, 2020 15:06

ayshih force-pushed the diffrotdynamic branch from f103660 to ff1cb6c Compare January 29, 2020 17:50

ayshih and others added 11 commits February 5, 2020 10:00

Override the size() property in SunPy frames to allow use as data-les…

678f0f7

…s coordinate frame attributes

Added RotatedSunFrame for defining coordinate frames that account for…

4d4d8c1

… solar rotation

Allow explicitly defining the coordinate frame for a WCS

ee77239

Addressed comments

1125e3e

Added an as_base() method to strip off the RotatedSunFrame layer

692575e

Added tests for RotatedSunFrame

d89c362

Added tutorial for RotatedSunFrame

42fb592

Adjusted doctests for more stable float comparisons and for readability

3f97a8f

Added upgrading of scalar base when duration is an array, and added r…

47d8108

…otated_time property

Added hypothesis-based testing

8169591

Addressed comments

bee98cd

ayshih force-pushed the diffrotdynamic branch from ff1cb6c to bee98cd Compare February 5, 2020 15:01

dstansby reviewed Feb 5, 2020

View reviewed changes

dstansby approved these changes Feb 5, 2020

View reviewed changes

Cadair merged commit 3c89110 into sunpy:master Feb 5, 2020

dstansby mentioned this pull request Feb 5, 2020

Improve differential rotation documentation #3766

Closed

3 tasks

ayshih mentioned this pull request Feb 10, 2020

Fixes to the transformation test for NorthOffsetFrame #3775

Merged

ayshih deleted the diffrotdynamic branch February 10, 2020 16:47

Cadair mentioned this pull request Mar 30, 2020

Reconcile new RotatedSunFrame with differential_rotate #3931

Closed

ayshih mentioned this pull request Apr 21, 2020

What's new for 2.0 #4040

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Coordinate frames that include solar differential rotation #3537

Coordinate frames that include solar differential rotation #3537

ayshih commented Nov 22, 2019 •

edited

Loading

pep8speaks commented Nov 22, 2019 •

edited

Loading

ghost commented Nov 22, 2019 •

edited by ghost

Loading

ayshih commented Nov 22, 2019

ayshih commented Nov 23, 2019

ayshih commented Nov 27, 2019 •

edited

Loading

ayshih commented Dec 18, 2019

ayshih commented Dec 18, 2019

Cadair commented Dec 19, 2019

Cadair commented Dec 19, 2019

ayshih commented Dec 19, 2019

ayshih commented Jan 17, 2020

Cadair commented Jan 17, 2020

Cadair left a comment •

edited

Loading

Cadair Jan 17, 2020

Cadair commented Jan 17, 2020

dstansby left a comment

dstansby Feb 5, 2020 •

edited

Loading

dstansby Feb 5, 2020

dstansby Feb 5, 2020

dstansby left a comment

Cadair commented Feb 5, 2020

dstansby commented Feb 5, 2020

		@@ -0,0 +1,390 @@
		.. _sunpy-coordinates-rotatedsunframe:

		Differential rotation using coordinate frames


		Creating coordinates
		^^^^^^^^^^^^^^^^^^^^

Coordinate frames that include solar differential rotation #3537

Coordinate frames that include solar differential rotation #3537

Conversation

ayshih commented Nov 22, 2019 • edited Loading

Walkthrough

Propagating a coordinate

Another example

Grid overlay on a Map

pep8speaks commented Nov 22, 2019 • edited Loading

Comment last updated at 2020-02-05 15:01:21 UTC

ghost commented Nov 22, 2019 • edited by ghost Loading

ayshih commented Nov 22, 2019

Multiple rotation durations

ayshih commented Nov 23, 2019

Transforming a coordinate to a frame with a different observation time

ayshih commented Nov 27, 2019 • edited Loading

Differentially rotating a map (versus the old approach)

ayshih commented Dec 18, 2019

ayshih commented Dec 18, 2019

Cadair commented Dec 19, 2019

Cadair commented Dec 19, 2019

ayshih commented Dec 19, 2019

ayshih commented Jan 17, 2020

Cadair commented Jan 17, 2020

Cadair left a comment • edited Loading

Choose a reason for hiding this comment

Cadair Jan 17, 2020

Choose a reason for hiding this comment

Cadair commented Jan 17, 2020

dstansby left a comment

Choose a reason for hiding this comment

dstansby Feb 5, 2020 • edited Loading

Choose a reason for hiding this comment

dstansby Feb 5, 2020

Choose a reason for hiding this comment

dstansby Feb 5, 2020

Choose a reason for hiding this comment

dstansby left a comment

Choose a reason for hiding this comment

Cadair commented Feb 5, 2020

dstansby commented Feb 5, 2020

ayshih commented Nov 22, 2019 •

edited

Loading

pep8speaks commented Nov 22, 2019 •

edited

Loading

ghost commented Nov 22, 2019 •

edited by ghost

Loading

ayshih commented Nov 27, 2019 •

edited

Loading

Cadair left a comment •

edited

Loading

dstansby Feb 5, 2020 •

edited

Loading