-
-
Notifications
You must be signed in to change notification settings - Fork 573
/
_transformations.py
1357 lines (1047 loc) · 54.3 KB
/
_transformations.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
Coordinate Transformation Functions
This module contains the functions for converting one
`sunpy.coordinates.frames` object to another.
.. warning::
The functions in this submodule should never be called directly, transforming
between coordinate frames should be done using the ``.transform_to`` methods
on `~astropy.coordinates.BaseCoordinateFrame` or
`~astropy.coordinates.SkyCoord` instances.
"""
import logging
from copy import deepcopy
from functools import wraps
import erfa
import numpy as np
import astropy.units as u
from astropy.constants import c as speed_of_light
from astropy.coordinates import (
HCRS,
ICRS,
ITRS,
BaseCoordinateFrame,
ConvertError,
HeliocentricMeanEcliptic,
get_body_barycentric,
)
from astropy.coordinates.baseframe import frame_transform_graph
from astropy.coordinates.builtin_frames import make_transform_graph_docs
from astropy.coordinates.builtin_frames.utils import get_jd12
from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
from astropy.coordinates.representation import (
CartesianRepresentation,
SphericalRepresentation,
UnitSphericalRepresentation,
)
from astropy.coordinates.transformations import FunctionTransform, FunctionTransformWithFiniteDifference
from astropy.time import Time
from sunpy import log
from sunpy.sun import constants
from sunpy.util.decorators import sunpycontextmanager
from .frames import (
_J2000,
GeocentricEarthEquatorial,
GeocentricSolarEcliptic,
GeocentricSolarMagnetospheric,
Geomagnetic,
Heliocentric,
HeliocentricEarthEcliptic,
HeliocentricInertial,
HeliographicCarrington,
HeliographicStonyhurst,
Helioprojective,
SolarMagnetic,
)
RSUN_METERS = constants.get('radius').si.to(u.m)
__all__ = ['transform_with_sun_center',
'propagate_with_solar_surface']
# Boolean flag for whether to ignore the motion of the center of the Sun in inertial space
_ignore_sun_motion = False
# If not None, the name of the differential-rotation model to use for any obstime change
_autoapply_diffrot = None
@sunpycontextmanager
def transform_with_sun_center():
"""
Context manager for coordinate transformations to ignore the motion of the center of the Sun.
Normally, coordinates refer to a point in inertial space (relative to the barycenter of the
solar system). Transforming to a different observation time does not move the point at all,
but rather only updates the coordinate representation as needed for the origin and axis
orientations at the new observation time. However, the center of the Sun moves over time.
Thus, for example, a coordinate that lies on the surface of the Sun at one observation time
will not continue to lie on the surface of the Sun at other observation times.
Under this context manager, transformations will instead move the coordinate over time to
"follow" the translational motion of the center of Sun, thus maintaining the position of the
coordinate relative to the center of the Sun.
Notes
-----
This context manager accounts only for the motion of the center of the Sun, i.e.,
translational motion. The motion of solar features due to any rotation of the Sun about its
rotational axis is not accounted for.
Due to the implementation approach, this context manager modifies transformations between only
these five coordinate frames:
`~sunpy.coordinates.frames.HeliographicStonyhurst`,
`~sunpy.coordinates.frames.HeliographicCarrington`,
`~sunpy.coordinates.frames.HeliocentricInertial`,
`~sunpy.coordinates.frames.Heliocentric`, and
`~sunpy.coordinates.frames.Helioprojective`.
Examples
--------
>>> from astropy.coordinates import SkyCoord
>>> from sunpy.coordinates import HeliographicStonyhurst, transform_with_sun_center
>>> import astropy.units as u
>>> start_frame = HeliographicStonyhurst(obstime="2001-01-01")
>>> end_frame = HeliographicStonyhurst(obstime="2001-02-01")
>>> sun_center = SkyCoord(0*u.deg, 0*u.deg, 0*u.AU, frame=start_frame)
>>> sun_center
<SkyCoord (HeliographicStonyhurst: obstime=2001-01-01T00:00:00.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, AU)
(0., 0., 0.)>
>>> sun_center.transform_to(end_frame) # transformations do not normally follow Sun center
<SkyCoord (HeliographicStonyhurst: obstime=2001-02-01T00:00:00.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, AU)
(23.33174233, -5.96399877, 0.00027959)>
>>> with transform_with_sun_center():
... sun_center.transform_to(end_frame) # now following Sun center
<SkyCoord (HeliographicStonyhurst: obstime=2001-02-01T00:00:00.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, AU)
(0., 0., 0.)>
"""
try:
global _ignore_sun_motion
old_ignore_sun_motion = _ignore_sun_motion # nominally False
if not old_ignore_sun_motion:
log.debug("Ignoring the motion of the center of the Sun for transformations")
_ignore_sun_motion = True
yield
finally:
if not old_ignore_sun_motion:
log.debug("Stop ignoring the motion of the center of the Sun for transformations")
_ignore_sun_motion = old_ignore_sun_motion
@sunpycontextmanager
def propagate_with_solar_surface(rotation_model='howard'):
"""
Context manager for coordinate transformations to automatically apply solar
differential rotation for any change in observation time.
Normally, coordinates refer to a point in inertial space (relative to the
barycenter of the solar system). Transforming to a different observation time
does not move the point at all, but rather only updates the coordinate
representation as needed for the origin and axis orientations at the new
observation time.
Under this context manager, transformations will instead treat the coordinate
as if it were referring to a point on the solar surface instead of a point in
inertial space. If a transformation has a change in observation time, the
heliographic longitude of the point will be updated according to the specified
rotation model.
Parameters
----------
rotation_model : `str`
Accepted model names are ``'howard'`` (default), ``'snodgrass'``,
``'allen'``, and ``'rigid'``. See the documentation for
:func:`~sunpy.sun.models.differential_rotation` for the differences
between these models.
Notes
-----
This context manager also ignores the motion of the center of the Sun (see
:func:`~sunpy.coordinates.transform_with_sun_center`).
Due to the implementation approach, this context manager modifies
transformations between only these five coordinate frames:
`~sunpy.coordinates.frames.HeliographicStonyhurst`,
`~sunpy.coordinates.frames.HeliographicCarrington`,
`~sunpy.coordinates.frames.HeliocentricInertial`,
`~sunpy.coordinates.frames.Heliocentric`, and
`~sunpy.coordinates.frames.Helioprojective`.
Examples
--------
.. minigallery:: sunpy.coordinates.propagate_with_solar_surface
>>> import astropy.units as u
>>> from astropy.coordinates import SkyCoord
>>> from sunpy.coordinates import HeliocentricInertial, propagate_with_solar_surface
>>> meridian = SkyCoord(0*u.deg, [-60, -30, 0, 30, 60]*u.deg, 1*u.AU,
... frame=HeliocentricInertial, obstime='2021-09-15')
>>> out_frame = HeliocentricInertial(obstime='2021-09-21')
>>> with propagate_with_solar_surface():
... print(meridian.transform_to(out_frame))
<SkyCoord (HeliocentricInertial: obstime=2021-09-21T00:00:00.000): (lon, lat, distance) in (deg, deg, AU)
[(70.24182965, -60., 1.),
(82.09298036, -30., 1.),
(85.9579703 , 0., 1.),
(82.09298036, 30., 1.),
(70.24182965, 60., 1.)]>
>>> with propagate_with_solar_surface(rotation_model='rigid'):
... print(meridian.transform_to(out_frame))
<SkyCoord (HeliocentricInertial: obstime=2021-09-21T00:00:00.000): (lon, lat, distance) in (deg, deg, AU)
[(85.1064, -60., 1.), (85.1064, -30., 1.),
(85.1064, 0., 1.), (85.1064, 30., 1.),
(85.1064, 60., 1.)]>
"""
with transform_with_sun_center():
try:
global _autoapply_diffrot
old_autoapply_diffrot = _autoapply_diffrot # nominally False
log.debug("Enabling automatic solar differential rotation "
f"('{rotation_model}') for any changes in obstime")
_autoapply_diffrot = rotation_model
yield
finally:
if not old_autoapply_diffrot:
log.debug("Disabling automatic solar differential rotation "
"for any changes in obstime")
_autoapply_diffrot = old_autoapply_diffrot
# Global counter to keep track of the layer of transformation
_layer_level = 0
def _transformation_debug(description):
"""
Decorator to produce debugging output for a transformation function: its description, inputs,
and output. Unicode box-drawing characters are used.
"""
def decorator(func):
@wraps(func)
def wrapped_func(*args, **kwargs):
global _layer_level
# Check if the logging level is at least DEBUG (for performance reasons)
debug_output = log.getEffectiveLevel() <= logging.DEBUG
if debug_output:
# Indentation for transformation layer
indentation = "\u2502 " * _layer_level
# For the input arguments, add indentation to any lines after the first line
from_str = str(args[0]).replace("\n", f"\n {indentation}\u2502 ")
to_str = str(args[1]).replace("\n", f"\n {indentation}\u2502 ")
# Log the description and the input arguments
log.debug(f"{indentation}{description}")
log.debug(f"{indentation}\u251c\u2500From: {from_str}")
log.debug(f"{indentation}\u251c\u2500To : {to_str}")
# Increment the layer level to increase the indentation for nested transformations
_layer_level += 1
result = func(*args, **kwargs)
if debug_output:
# Decrement the layer level
_layer_level -= 1
# For the output, add intention to any lines after the first line
out_str = str(result).replace("\n", f"\n {indentation} ")
# Log the output
log.debug(f"{indentation}\u2514\u2500Out : {out_str}")
return result
return wrapped_func
return decorator
def _observers_are_equal(obs_1, obs_2):
# Note that this also lets pass the situation where both observers are None
if obs_1 is obs_2:
return True
# obs_1 != obs_2
if obs_1 is None:
raise ConvertError("The source observer is set to None, but the transformation requires "
"the source observer to be specified, as the destination observer "
f"is set to {obs_2}.")
if obs_2 is None:
raise ConvertError("The destination observer is set to None, but the transformation "
"requires the destination observer to be specified, as the "
f"source observer is set to {obs_1}.")
if isinstance(obs_1, str):
if obs_1 == "self":
return False
raise ConvertError("The source observer needs to have `obstime` set because the "
"destination observer is different.")
if isinstance(obs_2, str):
if obs_2 == "self":
return False
raise ConvertError("The destination observer needs to have `obstime` set because the "
"source observer is different.")
return np.atleast_1d(u.allclose(obs_1.lat, obs_2.lat) and
u.allclose(obs_1.lon, obs_2.lon) and
u.allclose(obs_1.radius, obs_2.radius) and
_times_are_equal(obs_1.obstime, obs_2.obstime)).all()
def _check_observer_defined(frame):
if frame.observer is None:
raise ConvertError("This transformation cannot be performed because the "
f"{frame.__class__.__name__} frame has observer=None.")
elif isinstance(frame.observer, str):
if frame.observer != "self":
raise ConvertError("This transformation cannot be performed because the "
f"{frame.__class__.__name__} frame needs a specified obstime "
f"to fully resolve observer='{frame.observer}'.")
elif not isinstance(frame, HeliographicCarrington):
raise ConvertError(f"The {frame.__class__.__name__} frame has observer='self' "
"but this is valid for only HeliographicCarrington frames.")
def _times_are_equal(time_1, time_2):
# Checks whether times are equal
if isinstance(time_1, Time) and isinstance(time_2, Time):
# We explicitly perform the check in TAI to avoid possible numerical precision differences
# between a time in UTC and the same time after a UTC->TAI->UTC conversion
return np.all(time_1.tai == time_2.tai)
# We also deem the times equal if they are both None
return time_1 is None and time_2 is None
# =============================================================================
# ------------------------- Transformation Framework --------------------------
# =============================================================================
def _transform_obstime(frame, obstime):
"""
Transform a frame to a new obstime using the appropriate loopback transformation.
If the new obstime is None, no transformation is performed.
If the frame's obstime is None, the frame is copied with the new obstime.
"""
# If obstime is None or the obstime matches, nothing needs to be done
if obstime is None or _times_are_equal(frame.obstime, obstime):
return frame
# Transform to the new obstime using the appropriate loopback transformation
new_frame = frame.replicate(obstime=obstime)
if frame.obstime is not None:
return frame.transform_to(new_frame)
else:
return new_frame
def _rotation_matrix_hgs_to_hgc(obstime, observer_distance_from_sun):
"""
Return the rotation matrix from HGS to HGC at the same observation time
"""
if obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Import here to avoid a circular import
from .sun import L0, earth_distance
# Calculate the difference in light travel time if the observer is at a different distance from
# the Sun than the Earth is
delta_time = (observer_distance_from_sun - earth_distance(obstime)) / speed_of_light
# Calculate the corresponding difference in apparent longitude
delta_lon = delta_time * constants.sidereal_rotation_rate
# Rotation is only in longitude, so only around the Z axis
return rotation_matrix(-(L0(obstime) + delta_lon), 'z')
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicStonyhurst, HeliographicCarrington)
@_transformation_debug("HGS->HGC")
def hgs_to_hgc(hgscoord, hgcframe):
"""
Convert from Heliographic Stonyhurst to Heliographic Carrington.
"""
_check_observer_defined(hgcframe)
if isinstance(hgcframe.observer, str) and hgcframe.observer == "self":
observer_radius = hgscoord.radius
else:
observer_radius = hgcframe.observer.radius
# First transform the HGS coord to the HGC obstime
int_coord = _transform_obstime(hgscoord, hgcframe.obstime)
# Rotate from HGS to HGC
total_matrix = _rotation_matrix_hgs_to_hgc(int_coord.obstime, observer_radius)
newrepr = int_coord.cartesian.transform(total_matrix)
return hgcframe._replicate(newrepr, obstime=int_coord.obstime)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicCarrington, HeliographicStonyhurst)
@_transformation_debug("HGC->HGS")
def hgc_to_hgs(hgccoord, hgsframe):
"""
Convert from Heliographic Carrington to Heliographic Stonyhurst.
"""
_check_observer_defined(hgccoord)
hgccoord = hgccoord.make_3d()
if isinstance(hgccoord.observer, str) and hgccoord.observer == "self":
observer_radius = hgccoord.radius
else:
observer_radius = hgccoord.observer.radius
# First transform the HGC coord to the HGS obstime
int_coord = _transform_obstime(hgccoord, hgsframe.obstime)
# Rotate from HGC to HGS
total_matrix = matrix_transpose(_rotation_matrix_hgs_to_hgc(int_coord.obstime,
observer_radius))
newrepr = int_coord.cartesian.transform(total_matrix)
return hgsframe._replicate(newrepr, obstime=int_coord.obstime)
def _matrix_hcc_to_hpc():
# Returns the transformation matrix that permutes/swaps axes from HCC to HPC
# HPC spherical coordinates are a left-handed frame with these equivalent Cartesian axes:
# HPC_X = -HCC_Z
# HPC_Y = HCC_X
# HPC_Z = HCC_Y
# (HPC_X and HPC_Y are not to be confused with HPC_Tx and HPC_Ty)
return np.array([[0, 0, -1],
[1, 0, 0],
[0, 1, 0]])
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
Heliocentric, Helioprojective)
@_transformation_debug("HCC->HPC")
def hcc_to_hpc(helioccoord, heliopframe):
"""
Convert from Heliocentric Cartesian to Helioprojective Cartesian.
"""
_check_observer_defined(helioccoord)
_check_observer_defined(heliopframe)
# Transform the HPC observer (in HGS) to the HPC obstime in case it's different
observer = _transform_obstime(heliopframe.observer, heliopframe.obstime)
# Loopback transform HCC coord to obstime and observer of HPC frame
int_frame = Heliocentric(obstime=observer.obstime, observer=observer)
int_coord = helioccoord.transform_to(int_frame)
# Shift the origin from the Sun to the observer
distance = int_coord.observer.radius
newrepr = int_coord.cartesian - CartesianRepresentation(0*u.m, 0*u.m, distance)
# Permute/swap axes from HCC to HPC equivalent Cartesian
newrepr = newrepr.transform(_matrix_hcc_to_hpc())
# Explicitly represent as spherical because external code (e.g., wcsaxes) expects it
return heliopframe.realize_frame(newrepr.represent_as(SphericalRepresentation))
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
Helioprojective, Heliocentric)
@_transformation_debug("HPC->HCC")
def hpc_to_hcc(heliopcoord, heliocframe):
"""
Convert from Helioprojective Cartesian to Heliocentric Cartesian.
"""
_check_observer_defined(heliopcoord)
_check_observer_defined(heliocframe)
heliopcoord = heliopcoord.make_3d()
# Permute/swap axes from HPC equivalent Cartesian to HCC
newrepr = heliopcoord.cartesian.transform(matrix_transpose(_matrix_hcc_to_hpc()))
# Transform the HPC observer (in HGS) to the HPC obstime in case it's different
observer = _transform_obstime(heliopcoord.observer, heliopcoord.obstime)
# Shift the origin from the observer to the Sun
distance = observer.radius
newrepr += CartesianRepresentation(0*u.m, 0*u.m, distance)
# Complete the conversion of HPC to HCC at the obstime and observer of the HPC coord
int_coord = Heliocentric(newrepr, obstime=observer.obstime, observer=observer)
# Loopback transform HCC as needed
return int_coord.transform_to(heliocframe)
def _rotation_matrix_hcc_to_hgs(longitude, latitude):
# Returns the rotation matrix from HCC to HGS based on the observer longitude and latitude
# Permute the axes of HCC to match HGS Cartesian equivalent
# HGS_X = HCC_Z
# HGS_Y = HCC_X
# HGS_Z = HCC_Y
axes_matrix = np.array([[0, 0, 1],
[1, 0, 0],
[0, 1, 0]])
# Rotate in latitude and longitude (sign difference because of direction difference)
lat_matrix = rotation_matrix(latitude, 'y')
lon_matrix = rotation_matrix(-longitude, 'z')
return lon_matrix @ lat_matrix @ axes_matrix
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
Heliocentric, HeliographicStonyhurst)
@_transformation_debug("HCC->HGS")
def hcc_to_hgs(helioccoord, heliogframe):
"""
Convert from Heliocentric Cartesian to Heliographic Stonyhurst.
"""
_check_observer_defined(helioccoord)
# Transform the HCC observer (in HGS) to the HCC obstime in case it's different
hcc_observer_at_hcc_obstime = _transform_obstime(helioccoord.observer, helioccoord.obstime)
total_matrix = _rotation_matrix_hcc_to_hgs(hcc_observer_at_hcc_obstime.lon,
hcc_observer_at_hcc_obstime.lat)
# Transform from HCC to HGS at the HCC obstime
newrepr = helioccoord.cartesian.transform(total_matrix)
int_coord = HeliographicStonyhurst(newrepr, obstime=hcc_observer_at_hcc_obstime.obstime)
# For historical reasons, we support HCC with no obstime transforming to HGS with an obstime
if int_coord.obstime is None and heliogframe.obstime is not None:
int_coord = int_coord.replicate(obstime=heliogframe.obstime)
# Loopback transform HGS as needed
return int_coord.transform_to(heliogframe)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicStonyhurst, Heliocentric)
@_transformation_debug("HGS->HCC")
def hgs_to_hcc(heliogcoord, heliocframe):
"""
Convert from Heliographic Stonyhurst to Heliocentric Cartesian.
"""
_check_observer_defined(heliocframe)
heliogcoord = heliogcoord.make_3d()
# Loopback transform HGS if there is a change in obstime
int_coord = _transform_obstime(heliogcoord, heliocframe.obstime)
# Transform the HCC observer (in HGS) to the HCC obstime in case it's different
hcc_observer_at_hcc_obstime = _transform_obstime(heliocframe.observer, int_coord.obstime)
total_matrix = matrix_transpose(_rotation_matrix_hcc_to_hgs(hcc_observer_at_hcc_obstime.lon,
hcc_observer_at_hcc_obstime.lat))
# Transform from HGS to HCC at the same obstime
newrepr = int_coord.cartesian.transform(total_matrix)
return heliocframe.realize_frame(newrepr)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
Helioprojective, Helioprojective)
@_transformation_debug("HPC->HPC")
def hpc_to_hpc(from_coo, to_frame):
"""
This converts from HPC to HPC, with different observer location parameters.
It does this by transforming through HGS.
"""
if _observers_are_equal(from_coo.observer, to_frame.observer) and \
_times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
_check_observer_defined(from_coo)
_check_observer_defined(to_frame)
hgs = from_coo.transform_to(HeliographicStonyhurst(obstime=to_frame.obstime))
hpc = hgs.transform_to(to_frame)
return hpc
def _rotation_matrix_reprs_to_reprs(start_representation, end_representation):
"""
Return the matrix for the direct rotation from one representation to a second representation.
The representations need not be normalized first, and can be arrays of representations.
"""
A = start_representation.to_cartesian()
B = end_representation.to_cartesian()
rotation_axis = A.cross(B)
# Calculate the angle using both cross and dot products to minimize numerical-precision issues
rotation_angle = -np.arctan2(rotation_axis.norm(), A.dot(B)) # negation is required
if rotation_angle.isscalar:
# This line works around some input/output quirks of Astropy's rotation_matrix()
matrix = np.array(rotation_matrix(rotation_angle, rotation_axis.xyz.value.tolist()))
else:
matrix_list = [np.array(rotation_matrix(angle, axis.xyz.value.tolist()))
for angle, axis in zip(rotation_angle, rotation_axis)]
matrix = np.stack(matrix_list)
return matrix
def _rotation_matrix_reprs_to_xz_about_z(representations):
"""
Return one or more matrices for rotating one or more representations around the Z axis into the
XZ plane.
"""
A = representations.to_cartesian()
# Zero out the Z components
# (The additional transpose operations are to handle both scalar and array inputs)
A_no_z = CartesianRepresentation((A.xyz.T * [1, 1, 0]).T)
# Rotate the resulting vector to the X axis
x_axis = CartesianRepresentation(1, 0, 0)
matrix = _rotation_matrix_reprs_to_reprs(A_no_z, x_axis)
return matrix
def _sun_earth_icrf(time):
"""
Return the Sun-Earth vector for ICRF-based frames.
"""
sun_pos_icrs = get_body_barycentric('sun', time)
earth_pos_icrs = get_body_barycentric('earth', time)
return earth_pos_icrs - sun_pos_icrs
# The Sun's north pole is oriented RA=286.13 deg, dec=63.87 deg in ICRS, and thus HCRS as well
# (See Archinal et al. 2011,
# "Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009")
# The orientation of the north pole in ICRS/HCRS is assumed to be constant in time
_SOLAR_NORTH_POLE_HCRS = UnitSphericalRepresentation(lon=constants.get('alpha_0'),
lat=constants.get('delta_0'))
# Calculate the rotation matrix to de-tilt the Sun's rotation axis to be parallel to the Z axis
_SUN_DETILT_MATRIX = _rotation_matrix_reprs_to_reprs(_SOLAR_NORTH_POLE_HCRS,
CartesianRepresentation(0, 0, 1))
def _affine_params_hcrs_to_hgs(hcrs_time, hgs_time):
"""
Return the affine parameters (matrix and offset) from HCRS to HGS
HGS shares the same origin (the Sun) as HCRS, but has its Z axis aligned with the Sun's
rotation axis and its X axis aligned with the projection of the Sun-Earth vector onto the Sun's
equatorial plane (i.e., the component of the Sun-Earth vector perpendicular to the Z axis).
Thus, the transformation matrix is the product of the matrix to align the Z axis (by de-tilting
the Sun's rotation axis) and the matrix to align the X axis. The first matrix is independent
of time and is pre-computed, while the second matrix depends on the time-varying Sun-Earth
vector.
"""
# Determine the Sun-Earth vector in ICRS
# Since HCRS is ICRS with an origin shift, this is also the Sun-Earth vector in HCRS
sun_pos_icrs = get_body_barycentric('sun', hgs_time)
earth_pos_icrs = get_body_barycentric('earth', hgs_time)
sun_earth = earth_pos_icrs - sun_pos_icrs
# De-tilt the Sun-Earth vector to the frame with the Sun's rotation axis parallel to the Z axis
sun_earth_detilt = sun_earth.transform(_SUN_DETILT_MATRIX)
# Rotate the Sun-Earth vector about the Z axis so that it lies in the XZ plane
rot_matrix = _rotation_matrix_reprs_to_xz_about_z(sun_earth_detilt)
total_matrix = rot_matrix @ _SUN_DETILT_MATRIX
# All of the above is calculated for the HGS observation time
# If the HCRS observation time is different, calculate the translation in origin
if not _ignore_sun_motion and np.any(hcrs_time != hgs_time):
sun_pos_old_icrs = get_body_barycentric('sun', hcrs_time)
offset_icrf = sun_pos_old_icrs - sun_pos_icrs
else:
offset_icrf = sun_pos_icrs * 0 # preserves obstime shape
offset = offset_icrf.transform(total_matrix)
return total_matrix, offset
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HCRS, HeliographicStonyhurst)
@_transformation_debug("HCRS->HGS")
def hcrs_to_hgs(hcrscoord, hgsframe):
"""
Convert from HCRS to Heliographic Stonyhurst (HGS).
Even though we calculate the parameters for the affine transform, we use
``FunctionTransformWithFiniteDifference`` because otherwise there is no way to account for the
induced angular velocity when transforming a coordinate with velocity information.
"""
if hgsframe.obstime is None:
raise ConvertError("To perform this transformation, the HeliographicStonyhurst"
" frame needs a specified `obstime`.")
rot_matrix, offset = _affine_params_hcrs_to_hgs(hcrscoord.obstime, hgsframe.obstime)
return hgsframe.realize_frame(hcrscoord.cartesian.transform(rot_matrix) + offset)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicStonyhurst, HCRS)
@_transformation_debug("HGS->HCRS")
def hgs_to_hcrs(hgscoord, hcrsframe):
"""
Convert from Heliographic Stonyhurst to HCRS.
Even though we calculate the parameters for the affine transform, we use
``FunctionTransformWithFiniteDifference`` because otherwise there is no way to account for the
induced angular velocity when transforming a coordinate with velocity information.
"""
if hgscoord.obstime is None:
raise ConvertError("To perform this transformation, the HeliographicStonyhurst"
" frame needs a specified `obstime`.")
hgscoord = hgscoord.make_3d()
# Calculate the matrix and offset in the HCRS->HGS direction
forward_matrix, forward_offset = _affine_params_hcrs_to_hgs(hcrsframe.obstime, hgscoord.obstime)
# Invert the transformation to get the HGS->HCRS transformation
reverse_matrix = matrix_transpose(forward_matrix)
reverse_offset = (-forward_offset).transform(reverse_matrix)
return hcrsframe.realize_frame(hgscoord.cartesian.transform(reverse_matrix) + reverse_offset)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicStonyhurst, HeliographicStonyhurst)
@_transformation_debug("HGS->HGS")
def hgs_to_hgs(from_coo, to_frame):
"""
Convert between two Heliographic Stonyhurst frames.
"""
if to_frame.obstime is None:
return from_coo.replicate()
elif _times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
else:
if _autoapply_diffrot:
from_coo = from_coo._apply_diffrot((to_frame.obstime - from_coo.obstime).to('day'),
_autoapply_diffrot)
return from_coo.transform_to(HCRS(obstime=to_frame.obstime)).transform_to(to_frame)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicCarrington, HeliographicCarrington)
@_transformation_debug("HGC->HGC")
def hgc_to_hgc(from_coo, to_frame):
"""
Convert between two Heliographic Carrington frames.
"""
if _observers_are_equal(from_coo.observer, to_frame.observer) and \
_times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
_check_observer_defined(from_coo)
_check_observer_defined(to_frame)
# Convert through HGS
hgscoord = from_coo.transform_to(HeliographicStonyhurst(obstime=from_coo.obstime))
return hgscoord.transform_to(to_frame)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
Heliocentric, Heliocentric)
@_transformation_debug("HCC->HCC")
def hcc_to_hcc(from_coo, to_frame):
"""
Convert between two Heliocentric frames.
"""
if _observers_are_equal(from_coo.observer, to_frame.observer) and \
_times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
_check_observer_defined(from_coo)
_check_observer_defined(to_frame)
# Convert through HGS
hgscoord = from_coo.transform_to(HeliographicStonyhurst(obstime=to_frame.obstime))
return hgscoord.transform_to(to_frame)
def _rotation_matrix_hme_to_hee(hmeframe):
"""
Return the rotation matrix from HME to HEE at the same observation time
"""
# Get the Sun-Earth vector
sun_earth = HCRS(_sun_earth_icrf(hmeframe.obstime), obstime=hmeframe.obstime)
sun_earth_hme = sun_earth.transform_to(hmeframe).cartesian
# Rotate the Sun-Earth vector about the Z axis so that it lies in the XZ plane
rot_matrix = _rotation_matrix_reprs_to_xz_about_z(sun_earth_hme)
# Tilt the rotated Sun-Earth vector so that it is aligned with the X axis
tilt_matrix = _rotation_matrix_reprs_to_reprs(sun_earth_hme.transform(rot_matrix),
CartesianRepresentation(1, 0, 0))
return tilt_matrix @ rot_matrix
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliocentricMeanEcliptic, HeliocentricEarthEcliptic)
@_transformation_debug("HME->HEE")
def hme_to_hee(hmecoord, heeframe):
"""
Convert from Heliocentric Mean Ecliptic to Heliocentric Earth Ecliptic
"""
if heeframe.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Convert to the HME frame with mean equinox of date at the HEE obstime, through HCRS
int_frame = HeliocentricMeanEcliptic(obstime=heeframe.obstime, equinox=heeframe.obstime)
int_coord = hmecoord.transform_to(HCRS(obstime=hmecoord.obstime)).transform_to(int_frame)
# Rotate the intermediate coord to the HEE frame
total_matrix = _rotation_matrix_hme_to_hee(int_frame)
newrepr = int_coord.cartesian.transform(total_matrix)
return heeframe.realize_frame(newrepr)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliocentricEarthEcliptic, HeliocentricMeanEcliptic)
@_transformation_debug("HEE->HME")
def hee_to_hme(heecoord, hmeframe):
"""
Convert from Heliocentric Earth Ecliptic to Heliocentric Mean Ecliptic
"""
if heecoord.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
int_frame = HeliocentricMeanEcliptic(obstime=heecoord.obstime, equinox=heecoord.obstime)
# Rotate the HEE coord to the intermediate frame
total_matrix = matrix_transpose(_rotation_matrix_hme_to_hee(int_frame))
int_repr = heecoord.cartesian.transform(total_matrix)
int_coord = int_frame.realize_frame(int_repr)
# Convert to the HME frame through HCRS
return int_coord.transform_to(HCRS(obstime=int_coord.obstime)).transform_to(hmeframe)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliocentricEarthEcliptic, HeliocentricEarthEcliptic)
@_transformation_debug("HEE->HEE")
def hee_to_hee(from_coo, to_frame):
"""
Convert between two Heliocentric Earth Ecliptic frames.
"""
if _times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
elif to_frame.obstime is None:
return from_coo
else:
return from_coo.transform_to(HCRS(obstime=from_coo.obstime)).transform_to(to_frame)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliocentricEarthEcliptic, GeocentricSolarEcliptic)
@_transformation_debug("HEE->GSE")
def hee_to_gse(heecoord, gseframe):
"""
Convert from Heliocentric Earth Ecliptic to Geocentric Solar Ecliptic
"""
# First transform the HEE coord to the GSE obstime
int_coord = _transform_obstime(heecoord, gseframe.obstime)
if int_coord.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Import here to avoid a circular import
from .sun import earth_distance
# Find the Earth-object vector in the intermediate frame
sun_earth_int = earth_distance(int_coord.obstime) * CartesianRepresentation(1, 0, 0)
earth_object_int = int_coord.cartesian - sun_earth_int
# Flip the vector in X and Y, but leave Z untouched
# (The additional transpose operations are to handle both scalar and array inputs)
newrepr = CartesianRepresentation((earth_object_int.xyz.T * [-1, -1, 1]).T)
return gseframe._replicate(newrepr, obstime=int_coord.obstime)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
GeocentricSolarEcliptic, HeliocentricEarthEcliptic)
@_transformation_debug("GSE->HEE")
def gse_to_hee(gsecoord, heeframe):
"""
Convert from Geocentric Solar Ecliptic to Heliocentric Earth Ecliptic
"""
# First transform the GSE coord to the HEE obstime
int_coord = _transform_obstime(gsecoord, heeframe.obstime)
if int_coord.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Import here to avoid a circular import
from .sun import earth_distance
# Find the Sun-object vector in the intermediate frame
earth_sun_int = earth_distance(int_coord.obstime) * CartesianRepresentation(1, 0, 0)
sun_object_int = int_coord.cartesian - earth_sun_int
# Flip the vector in X and Y, but leave Z untouched
# (The additional transpose operations are to handle both scalar and array inputs)
newrepr = CartesianRepresentation((sun_object_int.xyz.T * [-1, -1, 1]).T)
return heeframe._replicate(newrepr, obstime=int_coord.obstime)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
GeocentricSolarEcliptic, GeocentricSolarEcliptic)
@_transformation_debug("GSE->GSE")
def gse_to_gse(from_coo, to_frame):
"""
Convert between two Geocentric Solar Ecliptic frames.
"""
if _times_are_equal(from_coo.obstime, to_frame.obstime):
return to_frame.realize_frame(from_coo.data)
else:
heecoord = from_coo.transform_to(HeliocentricEarthEcliptic(obstime=from_coo.obstime))
return heecoord.transform_to(to_frame)
def _rotation_matrix_hgs_to_hci(obstime):
"""
Return the rotation matrix from HGS to HCI at the same observation time
"""
z_axis = CartesianRepresentation(0, 0, 1)*u.m
if not obstime.isscalar:
z_axis = z_axis._apply('repeat', obstime.size)
# Get the ecliptic pole in HGS
ecliptic_pole = HeliocentricMeanEcliptic(z_axis, obstime=obstime, equinox=_J2000)
ecliptic_pole_hgs = ecliptic_pole.transform_to(HeliographicStonyhurst(obstime=obstime))
# Rotate the ecliptic pole to the -YZ plane, which aligns the solar ascending node with the X
# axis
rot_matrix = _rotation_matrix_reprs_to_xz_about_z(ecliptic_pole_hgs.cartesian)
xz_to_yz_matrix = rotation_matrix(-90*u.deg, 'z')
return xz_to_yz_matrix @ rot_matrix
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliographicStonyhurst, HeliocentricInertial)
@_transformation_debug("HGS->HCI")
def hgs_to_hci(hgscoord, hciframe):
"""
Convert from Heliographic Stonyhurst to Heliocentric Inertial
"""
hgscoord = hgscoord.make_3d()
# First transform the HGS coord to the HCI obstime
int_coord = _transform_obstime(hgscoord, hciframe.obstime)
if int_coord.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Rotate from HGS to HCI
total_matrix = _rotation_matrix_hgs_to_hci(int_coord.obstime)
newrepr = int_coord.cartesian.transform(total_matrix)
return hciframe._replicate(newrepr, obstime=int_coord.obstime)
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
HeliocentricInertial, HeliographicStonyhurst)
@_transformation_debug("HCI->HGS")
def hci_to_hgs(hcicoord, hgsframe):
"""
Convert from Heliocentric Inertial to Heliographic Stonyhurst
"""
# First transform the HCI coord to the HGS obstime
int_coord = _transform_obstime(hcicoord, hgsframe.obstime)
if int_coord.obstime is None:
raise ConvertError("To perform this transformation, the coordinate"
" frame needs a specified `obstime`.")
# Rotate from HCI to HGS
total_matrix = matrix_transpose(_rotation_matrix_hgs_to_hci(int_coord.obstime))
newrepr = int_coord.cartesian.transform(total_matrix)
return hgsframe._replicate(newrepr, obstime=int_coord.obstime)