Merge pull request #267 from dstansby/hmi-example

Hmi example

Author: dstansby

examples/using_pfsspy/plot_hmi.py

@@ -0,0 +1,79 @@
+"""
+HMI PFSS solutions
+------------------
+Calculating a PFSS solution from a HMI synoptic map.
+
+This example shows how to calcualte a PFSS solution from a HMI synoptic map.
+There are a couple of important things that this example shows:
+
+    - HMI maps have non-standard metadata, so this needs to be fixed
+    - HMI synoptic maps are very big (1440 x 3600), so need to be downsampled
+      in order to calculate the PFSS solution in a reasonable time.
+"""
+
+###############################################################################
+# First import the required modules
+import astropy.units as u
+import matplotlib.pyplot as plt
+from sunpy.net import Fido, attrs as a
+import sunpy.map
+import pfsspy
+import pfsspy.utils
+
+###############################################################################
+# Set up the search.
+#
+# Note that for SunPy versions earlier than 2.0, a time attribute is needed to
+# do the search, even if (in this case) it isn't used, as the synoptic maps are
+# labelled by Carrington rotation number instead of time
+time = a.Time('2010/01/01', '2010/01/01')
+series = a.jsoc.Series('hmi.synoptic_mr_polfil_720s')
+crot = a.jsoc.PrimeKey('CAR_ROT', 2210)
+
+###############################################################################
+# Do the search.
+#
+# If you use this code, please replace this email address
+# with your own one, registered here:
+# http://jsoc.stanford.edu/ajax/register_email.html
+result = Fido.search(time, series, crot, a.jsoc.Notify("[email protected]"))
+files = Fido.fetch(result)
+
+###############################################################################
+# Read in a file. This will read in the first file downloaded to a sunpy Map
+# object. Note that HMI maps have several bits of metadata that do not comply
+# to the FITS standard, so we need to fix them first.
+hmi_map = sunpy.map.Map(files[0])
+pfsspy.utils.fix_hmi_meta(hmi_map)
+print('Data shape: ', hmi_map.data.shape)
+
+###############################################################################
+# Since this map is far to big to calculate a PFSS solution quickly, lets
+# resample it down to a smaller size.
+hmi_map = hmi_map.resample([360, 180] * u.pix)
+print('New shape: ', hmi_map.data.shape)
+
+###############################################################################
+# Now calculate the PFSS solution
+nrho = 35
+rss = 2.5
+input = pfsspy.Input(hmi_map, nrho, rss)
+output = pfsspy.pfss(input)
+
+###############################################################################
+# Using the Output object we can plot the source surface field, and the
+# polarity inversion line.
+ss_br = output.source_surface_br
+# Create the figure and axes
+fig = plt.figure()
+ax = plt.subplot(projection=ss_br)
+
+# Plot the source surface map
+ss_br.plot()
+# Plot the polarity inversion line
+ax.plot_coord(output.source_surface_pils[0])
+# Plot formatting
+plt.colorbar()
+ax.set_title('Source surface magnetic field')
+
+plt.show()

pfsspy/utils.py

@@ -8,16 +8,20 @@
 
 import sunpy.io
 import sunpy.map
+import sunpy.time
+
+import pfsspy.map
 
 
 def fix_hmi_meta(hmi_map):
     """
     Fix non-compliant FITS metadata in HMI maps.
 
     This function:
-    - Corrects CUNIT2 from 'Sine Latitude' to 'deg'
-    - Corrects CDELT1 and CDELT2 for a CEA projection
-    - Populates the DATE-OBS keyword from the T_OBS keyword
+        - Corrects CUNIT2 from 'Sine Latitude' to 'deg'
+        - Corrects CDELT1 and CDELT2 for a CEA projection
+        - Populates the DATE-OBS keyword from the T_OBS keyword
+        - Sets the observer coordinate to the Earth
 
     Notes
     -----
@@ -43,7 +47,11 @@ def fix_hmi_meta(hmi_map):
         hmi_map.meta['cdelt1'] = np.abs(hmi_map.meta['cdelt1'])
 
     if 'date-obs' not in hmi_map.meta and 't_obs' in hmi_map.meta:
-        hmi_map.meta['date-obs'] = hmi_map.meta['t_obs']
+        hmi_map.meta['date-obs'] = sunpy.time.parse_time(hmi_map.meta['t_obs']).isot
+
+    # Fix observer coordinate
+    if 'hglt_obs' not in hmi_map.meta:
+        hmi_map.meta.update(pfsspy.map._earth_obs_coord_meta(hmi_map.meta['date-obs']))
 
 
 def load_adapt(adapt_path):