Merge pull request #206 from dstansby/ex-improve

0.5.3 backports

doc/source/changes.rst

@@ -1,6 +1,14 @@
 Changelog
 =========
 
+0.5.3
+-----
+- Improved descriptions in the AIA overplotting example.
+- Fixed the 'date-obs' keyword in GONG metadata. Previously this just stored
+  the date and not the time; now both the date and time are properly stored.
+- Drastically sped up the calculation of source surface and solar surface
+  magnetic field footpoints.
+
 0.5.2
 -----
 - Fixed a bug in the GONG synoptic map source where a map failed to load once

examples/finding_data/adapt_MapSequence.py

@@ -2,25 +2,28 @@
 Parsing ADAPT Ensemble .fits files
 ==================================
 
-Parse an ADAPT FITS file into a `sunpy.map.MapSequence`. 
+Parse an ADAPT FITS file into a `sunpy.map.MapSequence`.
 """
 
 ###############################################################################
 # Necessary imports
 from adapt_helpers import example_adapt_map
-import sunpy.map,sunpy.io
-import matplotlib.pyplot as plt, matplotlib.gridspec as gridspec
+import sunpy.map
+import sunpy.io
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
 
 ###############################################################################
 # Load an example ADAPT fits file, utility stored in adapt_helpers.py
 adapt_fname = example_adapt_map()
 
 ###############################################################################
-# ADAPT synoptic magnetograms contain 12 realizations of synoptic magnetograms 
-# output as a result of varying model assumptions. See [here](https://www.swpc.noaa.gov/sites/default/files/images/u33/SWW_2012_Talk_04_27_2012_Arge.pdf))
-# 
-# Because the fits data is 3D, it cannot be passed directly to `sunpy.map.Map`, 
-# because this will take the first slice only and the other realizations are 
+# ADAPT synoptic magnetograms contain 12 realizations of synoptic magnetograms
+# output as a result of varying model assumptions. See
+# [here](https://www.swpc.noaa.gov/sites/default/files/images/u33/SWW_2012_Talk_04_27_2012_Arge.pdf))
+#
+# Because the fits data is 3D, it cannot be passed directly to `sunpy.map.Map`,
+# because this will take the first slice only and the other realizations are
 # lost. We want to end up with a `sunpy.map.MapSequence` containing all these
 # realiations as individual maps. These maps can then be individually accessed
 # and PFSS solutions generated from them.
@@ -29,29 +32,29 @@
 adapt_fits = sunpy.io.fits.read(adapt_fname)
 
 ###############################################################################
-# `adapt_fits` is a list of `HDPair` objects. The first of these contains the 
+# `adapt_fits` is a list of `HDPair` objects. The first of these contains the
 # 12 realizations data and a header with sufficient information to build the
 # `MapSequence`. We unpack this `HDPair` into a list of `(data,header)` tuples
 # where `data` are the different adapt realizations.
-data_header_pairs = [(map_slice,adapt_fits[0].header) 
-					 for map_slice in adapt_fits[0].data
-					 ]
+data_header_pairs = [(map_slice, adapt_fits[0].header)
+				     for map_slice in adapt_fits[0].data]
+
 
 ###############################################################################
 # Next, pass this list of tuples as the argument to sunpy.map.Map to create the
 # map sequence :
-adaptMapSequence = sunpy.map.Map(data_header_pairs,sequence=True) 
+adaptMapSequence = sunpy.map.Map(data_header_pairs, sequence=True)
 
 ###############################################################################
 # `adapt_map_sequence` is now a list of our individual adapt realizations. Note
 # the `.peek()` and `.plot()` methods of `MapSequence` returns instances of
-# `sunpy.visualization.MapSequenceAnimator` and 
+# `sunpy.visualization.MapSequenceAnimator` and
 # `matplotlib.animation.FuncAnimation1` instances. Here, we generate a static
 # plot accessing the individual maps in turn :
-fig = plt.figure(figsize=(7,8)) 
-gs = gridspec.GridSpec(4,3,figure=fig)
-for ii,aMap in enumerate(adaptMapSequence) : 
-    ax=fig.add_subplot(gs[ii],projection=aMap) 
-    aMap.plot(axes=ax,cmap='bwr',vmin=-2,vmax=2,title=f"Realization {1+ii:02d}") 
-plt.tight_layout(pad=5,h_pad=2)
-plt.show()   
+fig = plt.figure(figsize=(7, 8))
+gs = gridspec.GridSpec(4, 3, figure=fig)
+for ii, aMap in enumerate(adaptMapSequence):
+	ax = fig.add_subplot(gs[ii], projection=aMap)
+	aMap.plot(axes=ax, cmap='bwr', vmin=-2, vmax=2, title=f"Realization {1+ii:02d}")
+plt.tight_layout(pad=5, h_pad=2)
+plt.show()

examples/using_pfsspy/plot_aia_overplotting.py

@@ -82,12 +82,22 @@
 
 
 ###############################################################################
-# Now we construct a 10 x 10 grid of footpoitns to trace some magnetic field
-# lines from.
-s, phi = np.meshgrid(np.linspace(0.1, 0.2, 5),
-                     np.deg2rad(np.linspace(55, 65, 5)))
+# Now we construct a 5 x 5 grid of footpoitns to trace some magnetic field
+# lines from. These coordinates are defined in the native Carrington
+# coordinates of the input magnetogram.
+
+# Create 5 points spaced between sin(lat)={0.1, 0.2}
+s = np.linspace(0.1, 0.2, 5)
+# Create 5 points spaced between long={55, 65} degrees
+phi = np.linspace(55, 65, 5)
+print(f's = {s}')
+print(f'phi = {phi}')
+# Make a 2D grid from these 1D points
+s, phi = np.meshgrid(s, phi)
+
+# Now convert the points to a coordinate object
 lat = np.arcsin(s) * u.rad
-lon = phi * u.rad
+lon = phi * u.deg
 seeds = SkyCoord(lon.ravel(), lat.ravel(), 1.01 * const.R_sun,
                  frame=gong_map.coordinate_frame)
 
@@ -102,6 +112,9 @@
 
 ax.plot_coord(seeds, color='black', marker='o', linewidth=0, markersize=2)
 
+# Set the axes limits. These limits have to be in pixel values
+ax.set_xlim(0, 180)
+ax.set_ylim(45, 135)
 ax.set_title('Field line footpoints')
 ax.set_ylim(bottom=0)
 
@@ -126,16 +139,16 @@
 for fline in flines:
     ax.plot_coord(fline.coords, color='black', linewidth=1)
 
-# ax.set_xlim(55, 65)
-# ax.set_ylim(0.1, 0.25)
+# Set the axes limits. These limits have to be in pixel values
+ax.set_xlim(0, 180)
+ax.set_ylim(45, 135)
 ax.set_title('Photospheric field and traced field lines')
 ###############################################################################
 # Plot the AIA map, along with the traced magnetic field lines. Inside the
 # loop the field lines are converted to the AIA observer coordinate frame,
 # and then plotted on top of the map.
 fig = plt.figure()
 ax = plt.subplot(1, 1, 1, projection=aia)
-transform = ax.get_transform('world')
 aia.plot(ax)
 for fline in flines:
     ax.plot_coord(fline.coords, alpha=0.8, linewidth=1, color='black')

pfsspy/fieldline.py

@@ -84,30 +84,23 @@ def source_surface_feet(self):
         """
         Coordinates of the source suface footpoints.
         """
-        source_surface_feet = [fline.source_surface_footpoint for
-                               fline in self.field_lines]
+        x = np.array([fline._x[fline._ss_coord_index] for fline in self.field_lines])
+        y = np.array([fline._y[fline._ss_coord_index] for fline in self.field_lines])
+        z = np.array([fline._z[fline._ss_coord_index] for fline in self.field_lines])
 
-        if len(source_surface_feet) == 1:
-            source_surface_feet = source_surface_feet[0]
-        else:
-            source_surface_feet = coord.concatenate(source_surface_feet)
-
-        return source_surface_feet
+        return FieldLine._coords(x, y, z, self.field_lines[0]._output)
 
     @property
     @functools.lru_cache()
     def solar_feet(self):
         """
         Coordinates of the solar footpoints.
         """
-        solar_feet = [fline.solar_footpoint for fline in self.field_lines]
-
-        if len(solar_feet) == 1:
-            solar_feet = solar_feet[0]
-        else:
-            solar_feet = coord.concatenate(solar_feet)
+        x = np.array([fline._x[fline._solar_coord_index] for fline in self.field_lines])
+        y = np.array([fline._y[fline._solar_coord_index] for fline in self.field_lines])
+        z = np.array([fline._z[fline._solar_coord_index] for fline in self.field_lines])
 
-        return solar_feet
+        return FieldLine._coords(x, y, z, self.field_lines[0]._output)
 
 
 class ClosedFieldLines(FieldLines):
@@ -148,12 +141,15 @@ def coords(self):
         """
         Field line `~astropy.coordinates.SkyCoord`.
         """
-        r, lat, lon = coord.cartesian_to_spherical(
-            self._x, self._y, self._z)
+        return self._coords(self._x, self._y, self._z, self._output)
+
+    @staticmethod
+    def _coords(x, y, z, output):
+        r, lat, lon = coord.cartesian_to_spherical(x, y, z)
         r *= const.R_sun
-        lon += self._output._lon0 + 180 * u.deg
+        lon += output._lon0 + 180 * u.deg
         coords = coord.SkyCoord(
-            lon, lat, r, frame=self._output.coordinate_frame)
+            lon, lat, r, frame=output.coordinate_frame)
         return coords
 
     @property
@@ -195,10 +191,7 @@ def solar_footpoint(self):
         method returns the field line pointing out from the solar surface in
         this case.
         """
-        if self.polarity == 1 or not self.is_open:
-            return self.coords[0]
-        else:
-            return self.coords[-1]
+        return self.coords[self._solar_coord_index]
 
     @property
     def source_surface_footpoint(self):
@@ -218,10 +211,22 @@ def source_surface_footpoint(self):
         method returns the field line pointing out from the solar surface in
         this case.
         """
+        return self.coords[self._ss_coord_index]
+
+    @property
+    def _ss_coord_index(self):
+        """
+        Return 0 or -1 depending on which end of the coordinate array is the
+        source surface footpoint.
+        """
         if self.polarity == 1 or not self.is_open:
-            return self.coords[-1]
+            return -1
         else:
-            return self.coords[0]
+            return 0
+
+    @property
+    def _solar_coord_index(self):
+        return -1 - self._ss_coord_index
 
     @property
     @functools.lru_cache()
@@ -242,12 +247,8 @@ def expansion_factor(self):
         if not self.is_open:
             return np.nan
 
-        if self._r[0] > self._r[-1]:
-            solar_foot = -1
-            source_foot = 0
-        else:
-            solar_foot = 0
-            source_foot = -1
+        solar_foot = self._solar_coord_index
+        source_foot = self._ss_coord_index
 
         def interp(map, idx):
             x, y, z = self._x[idx], self._y[idx], self._z[idx]

pfsspy/map.py

@@ -16,6 +16,9 @@ def __init__(self, data, header, **kwargs):
                 # that value (see Thompson 2005)
                 header['CDELT2'] = 180 / np.pi * header['CDELT2']
                 header['KEYCOMMENTS']['CDELT2'] = '180 / pi * sine-lat step'
+            if 'time-obs' in header:
+                header['date-obs'] = (header['date-obs'] + ' ' +
+                                      header.pop('time-obs'))
         super().__init__(data, header, **kwargs)
 
     @classmethod

pfsspy/tests/test_map.py

@@ -28,6 +28,7 @@ def test_gong_source(gong_map):
 
     # Check round-trip is robust against sunpy changes to the meta
     m = sunpy.map.Map(m.data, m.meta)
+    assert m.date.isot == '2019-03-10T00:14:00.000'
 
 
 @pytest.fixture

pfsspy/tests/test_pfss.py

@@ -144,8 +144,9 @@ def test_header_generation():
     np.testing.assert_almost_equal(
         header['CDELT2'], (180 / np.pi) * (2 / ntheta))
 
-    assert header['CRPIX1'] == (nphi / 2) + 0.5
-    assert header['CRPIX2'] == (ntheta / 2) + 0.5
+    # Extra + 1s are for FITS counting from 1 indexing
+    assert header['CRPIX1'] == (nphi / 2) + 0.5 + 1
+    assert header['CRPIX2'] == (ntheta / 2) + 0.5 + 1
     assert header['CRVAL1'] == 0
     assert header['CRVAL2'] == 0
     assert header['CUNIT1'] == 'deg'

requirements/docs.txt

@@ -9,3 +9,4 @@ sphinx>2
 bs4
 drms
 zeep
+importlib_metadata