diff --git a/doc/source/changes.rst b/doc/source/changes.rst index bafd5159..5e77008a 100644 --- a/doc/source/changes.rst +++ b/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 diff --git a/examples/finding_data/adapt_MapSequence.py b/examples/finding_data/adapt_MapSequence.py index b0fcc293..3a7e5603 100644 --- a/examples/finding_data/adapt_MapSequence.py +++ b/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() diff --git a/examples/using_pfsspy/plot_aia_overplotting.py b/examples/using_pfsspy/plot_aia_overplotting.py index 6fadb3d1..a2fa534b 100644 --- a/examples/using_pfsspy/plot_aia_overplotting.py +++ b/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,8 +139,9 @@ 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 @@ -135,7 +149,6 @@ # 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') diff --git a/pfsspy/fieldline.py b/pfsspy/fieldline.py index e5e548cd..a3293e74 100644 --- a/pfsspy/fieldline.py +++ b/pfsspy/fieldline.py @@ -84,15 +84,11 @@ 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() @@ -100,14 +96,11 @@ 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] diff --git a/pfsspy/map.py b/pfsspy/map.py index 759840c1..964e0aba 100644 --- a/pfsspy/map.py +++ b/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 diff --git a/pfsspy/tests/test_map.py b/pfsspy/tests/test_map.py index d45754af..39b57937 100644 --- a/pfsspy/tests/test_map.py +++ b/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 diff --git a/pfsspy/tests/test_pfss.py b/pfsspy/tests/test_pfss.py index a2fa640b..5bec8bc2 100644 --- a/pfsspy/tests/test_pfss.py +++ b/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' diff --git a/requirements/docs.txt b/requirements/docs.txt index 16381102..e8edb9eb 100644 --- a/requirements/docs.txt +++ b/requirements/docs.txt @@ -9,3 +9,4 @@ sphinx>2 bs4 drms zeep +importlib_metadata