Merge pull request #332 from dstansby/open-closed

Plot open/closed map in the examples

doc/source/changes.rst

@@ -3,12 +3,21 @@
 Changelog
 =========
 
+1.0.1
+-----
+Bug fixes
+~~~~~~~~~
+- Fixed compatibility of map validity checks with sunpy 3.1.
+- Updated this changelog to make it clear that pfsspy 1.0.0 depends on
+  sunpy >= 3.0.
+
 1.0.0
 -----
 
 New requirements
 ~~~~~~~~~~~~~~~~
-pfsspy now depends on python >= 3.7, and now does *not* depend on Matplotlib.
+pfsspy now depends on python >= 3.7, sunpy >=3,
+and now does *not* depend on Matplotlib.
 
 New features
 ~~~~~~~~~~~~

examples/using_pfsspy/open_closed_map.py → examples/using_pfsspy/plot_open_closed_map.py

@@ -26,7 +26,7 @@
 
 ###############################################################################
 # Set the model parameters
-nrho = 60
+nrho = 40
 rss = 2.5
 
 ###############################################################################
@@ -44,7 +44,7 @@
 
 r = const.R_sun
 # Number of steps in cos(latitude)
-nsteps = 90
+nsteps = 45
 lon_1d = np.linspace(0, 2 * np.pi, nsteps * 2 + 1)
 lat_1d = np.arcsin(np.linspace(-1, 1, nsteps + 1))
 lon, lat = np.meshgrid(lon_1d, lat_1d, indexing='ij')

pfsspy/fieldline.py

@@ -130,8 +130,12 @@ def __init__(self, x, y, z, output):
         # Field line is open if one end is on the solar surface and one on
         # the source surface
         atol = 0.1
-        self._is_open = np.abs(self._r[0] - self._r[-1]) > atol
-        self._polarity = -np.sign(self._r[0] - self._r[-1]) * self._is_open
+        if len(self._r) <= 1:
+            self._is_open = False
+            self._polarity = 0
+        else:
+            self._is_open = np.abs(self._r[0] - self._r[-1]) > atol
+            self._polarity = -np.sign(self._r[0] - self._r[-1]) * self._is_open
 
     def __len__(self):
         return len(self._x)

pfsspy/utils.py

@@ -195,18 +195,18 @@ def is_full_sun_synoptic_map(m, error=False):
 def _is_full_sun_car(m, error=False):
     shape = m.data.shape
 
-    dphi = m.meta['cdelt1']
-    phi = shape[1] * dphi
-    if not np.allclose(phi, 360, atol=0.1):
+    dphi = m.scale.axis1
+    phi = shape[1] * u.pix * dphi
+    if not np.allclose(np.abs(phi), 360 * u.deg, atol=0.1 * u.deg):
         if error:
             raise ValueError('Number of points in phi direction times '
                              'CDELT1 must be close to 360 degrees. '
                              f'Instead got {dphi} x {shape[0]} = {phi}')
         return False
 
-    dtheta = m.meta['cdelt2']
-    theta = shape[0] * dtheta
-    if not np.allclose(theta, 180, atol=0.1):
+    dtheta = m.scale.axis2
+    theta = shape[0] * u.pix * dtheta
+    if not np.allclose(theta, 180 * u.deg, atol=0.1 * u.deg):
         if error:
             raise ValueError('Number of points in theta direction times '
                              'CDELT2 must be close to 180 degrees. '
@@ -218,18 +218,18 @@ def _is_full_sun_car(m, error=False):
 def _is_full_sun_cea(m, error=False):
     shape = m.data.shape
 
-    dphi = m.meta['cdelt1']
-    phi = shape[1] * dphi
-    if not np.allclose(phi, 360, atol=0.1):
+    dphi = m.scale.axis1
+    phi = shape[1] * u.pix * dphi
+    if not np.allclose(np.abs(phi), 360 * u.deg, atol=0.1 * u.deg):
         if error:
             raise ValueError('Number of points in phi direction times '
                              'CDELT1 must be close to 360 degrees. '
                              f'Instead got {dphi} x {shape[1]} = {phi}')
         return False
 
-    dtheta = m.meta['cdelt2']
-    theta = shape[0] * dtheta * np.pi / 2
-    if not np.allclose(theta, 180, atol=0.1):
+    dtheta = m.scale.axis2
+    theta = shape[0] * u.pix * dtheta * np.pi / 2
+    if not np.allclose(theta, 180 * u.deg, atol=0.1 * u.deg):
         if error:
             raise ValueError('Number of points in theta direction times '
                              'CDELT2 times pi/2 must be close to '

setup.cfg

@@ -36,7 +36,7 @@ install_requires =
   numpy
   scikit-image
   scipy
-  sunpy>=3
+  sunpy>=3,!=3.1.0
 
 [options.extras_require]
 docs =