Merge branch 'main' into solar_rotate_fix

CHANGELOG.rst

@@ -1060,7 +1060,7 @@ Bug Fixes
 Added/Improved Documentation
 ----------------------------
 
-- Added a gallery example (:ref:`sphx_glr_generated_gallery_plotting_plot_rectangle.py`) for drawing rectangles on maps. (`#4528 <https://github.com/sunpy/sunpy/pull/4528>`__)
+- Added a gallery example for drawing rectangles on maps. (`#4528 <https://github.com/sunpy/sunpy/pull/4528>`__)
 - Added an example (:ref:`sphx_glr_generated_gallery_plotting_wcsaxes_plotting_example.py`)
   of how pixel and SkyCoords work when plotted with `~astropy.visualization.wcsaxes`. (`#4867 <https://github.com/sunpy/sunpy/pull/4867>`__)
 - Added a gallery example  (:ref:`sphx_glr_generated_gallery_plotting_plotting_blank_map.py`) on how to create a blank map and mark locations. (`#5077 <https://github.com/sunpy/sunpy/pull/5077>`__)

changelog/6736.feature.rst

@@ -0,0 +1 @@
+:func:`sunpy.io.read_file` will now try to detect the filetype based on the content and then fallback to using the file extension.

changelog/7432.feature.rst

@@ -0,0 +1 @@
+Update ASDF schemas for upcoming ASDF standard 1.6.0.

changelog/7435.doc.rst

@@ -0,0 +1 @@
+Convert draw rectangle gallery example into a how-to guide(:ref:`sunpy-how-to-create-rectangle-on-map`)

changelog/7530.bugfix.rst

@@ -0,0 +1 @@
+Fixed an inaccuracy in the implementation of `~sunpy.coordinates.HeliocentricEarthEcliptic` and `~sunpy.coordinates.GeocentricSolarEcliptic` such that the Earth was not exactly in the XY plane, but rather had an error of up ~10 meters.

docs/how_to/create_rectangle_on_map.rst

@@ -0,0 +1,97 @@
+.. _sunpy-how-to-create-rectangle-on-map:
+
+*********************************************
+How to draw a rectangle on a `~sunpy.map.Map`
+*********************************************
+
+``sunpy`` provides a convenient method called :meth:`~sunpy.map.GenericMap.draw_quadrangle` to draw rectangles on maps.
+In this guide, we will demonstrate four different methods to draw a rectangle on a `sunpy.map.Map`.
+
+Specify corners with a single `~astropy.coordinates.SkyCoord`
+=============================================================
+
+We will use one `~astropy.coordinates.SkyCoord` to represent the two opposite corners.
+
+.. plot::
+    :include-source:
+    :context: close-figs
+
+    import matplotlib.pyplot as plt
+
+    import astropy.units as u
+    from astropy.coordinates import SkyCoord
+
+    import sunpy.data.sample
+    import sunpy.map
+
+    aia_map = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)
+
+    fig = plt.figure()
+    ax = fig.add_subplot(projection=aia_map)
+
+    aia_map.plot(axes=ax, clip_interval=(1, 99.99)*u.percent)
+
+    coords = SkyCoord(Tx=(100, 500) * u.arcsec, Ty=(200, 500) * u.arcsec,frame=aia_map.coordinate_frame)
+    aia_map.draw_quadrangle(coords, axes=ax, edgecolor="blue")
+
+    plt.show()
+
+Specify corners with separate `~astropy.coordinates.SkyCoord`
+=============================================================
+
+We will use two `~astropy.coordinates.SkyCoord` to represent the two opposite corners.
+
+.. plot::
+    :include-source:
+    :context: close-figs
+
+    fig = plt.figure()
+    ax = fig.add_subplot(projection=aia_map)
+
+    aia_map.plot(axes=ax, clip_interval=(1, 99.99)*u.percent)
+
+    bottom_left = SkyCoord(100 * u.arcsec, 200 * u.arcsec, frame=aia_map.coordinate_frame)
+    top_right = SkyCoord(500 * u.arcsec, 500 * u.arcsec, frame=aia_map.coordinate_frame)
+    aia_map.draw_quadrangle(bottom_left, axes=ax, top_right=top_right, edgecolor="green")
+
+    plt.show()
+
+Specify one corner with a width and height
+==========================================
+
+We will use one `~astropy.coordinates.SkyCoord` to represent the bottom left and supply a width and height to complete the rectangle.
+
+.. plot::
+    :include-source:
+    :context: close-figs
+
+    fig = plt.figure()
+    ax = fig.add_subplot(projection=aia_map)
+
+    aia_map.plot(axes=ax, clip_interval=(1, 99.99)*u.percent)
+
+    bottom_left = SkyCoord(100 * u.arcsec, 200 * u.arcsec, frame=aia_map.coordinate_frame)
+    width = 400 * u.arcsec
+    height = 300 * u.arcsec
+    aia_map.draw_quadrangle(bottom_left, axes=ax, width=width, height=height, edgecolor="yellow")
+
+    plt.show()
+
+Using pixel coordinates
+=======================
+
+We will use a `~astropy.coordinates.SkyCoord` to work out the pixel coordinates instead of using coordinates as we do above.
+
+.. plot::
+    :include-source:
+    :context: close-figs
+
+    fig = plt.figure()
+    ax = fig.add_subplot(projection=aia_map)
+    aia_map.plot(axes=ax, clip_interval=(1, 99.99)*u.percent)
+
+    bottom_left = aia_map.wcs.pixel_to_world(551 * u.pixel, 594 * u.pixel)
+    top_right = aia_map.wcs.pixel_to_world(717 * u.pixel, 719 * u.pixel)
+    aia_map.draw_quadrangle(bottom_left, axes=ax, top_right=top_right, edgecolor="red")
+
+    plt.show()

docs/how_to/index.rst

@@ -13,6 +13,7 @@ If you're starting fresh you might want to check out the :ref:`sunpy-tutorial-in
 
    coord_components
    create_a_map
+   create_rectangle_on_map
    create_coords
    create_custom_map
    create_custom_timeseries

sunpy/coordinates/_transformations.py

@@ -589,7 +589,8 @@ def _rotation_matrix_reprs_to_reprs(start_representation, end_representation):
     A = start_representation.to_cartesian()
     B = end_representation.to_cartesian()
     rotation_axis = A.cross(B)
-    rotation_angle = -np.arccos(A.dot(B) / (A.norm() * B.norm()))  # negation is required
+    # 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()

sunpy/coordinates/tests/test_transformations.py

@@ -695,6 +695,13 @@ def test_hme_hee_sunspice():
     assert_quantity_allclose(new.distance, old.distance)
 
 
+def test_hee_earth():
+    # The Earth in HEE should have negligible Z component
+    times = parse_time('2013-08-10 12:00') + np.arange(10) * u.s
+    earth_hee = get_earth(times).heliocentricearthecliptic
+    assert_quantity_allclose(0*u.m, earth_hee.cartesian.z, atol=1e-4*u.m)
+
+
 def test_hee_hee():
     # Test HEE loopback transformation
     obstime = Time('2001-01-01')