Solar Magnetic Field Extrapolation¶

This is an affiliated package for the AstroPy package. The documentation for this package is here:

SunPy Solar Magnetic Field Extrapolation Gallery¶

It does what it says on the tin.

Generating Example Gaussian Boundary Data

Generating Example Gaussian Boundary Data¶

In this example you will be generating some example data and extrapolate this using the basic potential extrapolator.

You can start by importing the necessary module components.

# Module imports
from solarbextrapolation.example_data_generator import generate_example_data, dummyDataToMap

You also need the ability to convert astropyunits.

import astropy.units as u

You need to define the parameters of the eare, includsing the x and y ranges as astropy quantities with angular or distance units and the grid shape.

# Input parameters:
arr_grid_shape = [ 20, 22 ]         # [ y-size, x-size ]
qua_xrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
qua_yrange = u.Quantity([ -11.0, 11.0 ] * u.arcsec)

# To randomly generate 2 poles simply don't add any pole parameters:
arr_Data = generate_example_data(arr_grid_shape, qua_xrange, qua_yrange)
# Note: each time you run this pole positions/magnitudes will change.

We can now convert this into a a sunpy map object:

aMap = dummyDataToMap(arr_Data, qua_xrange, qua_yrange)

We can see this map using peek:

aMap.peek()

arrA0 = [ Position, size, Max Magnitude ]

arrA0 = [ u.Quantity([ 25, 25 ] * u.percent), 10.0 * u.percent,  0.2 * u.T ]
arrA1 = [ u.Quantity([ 75, 75 ] * u.percent), 10.0 * u.percent, -0.2 * u.T ]

# To generate and view:
arr_Data = generate_example_data(arr_grid_shape, qua_xrange, qua_yrange, arrA0, arrA1)
aMap = dummyDataToMap(arr_Data, qua_xrange, qua_yrange)
aMap.peek()

auto_examplesimagessphx_glr_plot_gaussian_example_data_002.png

But absolute positioning using the map range units is also possible

arrA2 = [ u.Quantity([ -6,  6 ] * u.arcsec), 2 * u.arcsec, -0.2 * u.T ]
arrA3 = [ u.Quantity([  6, -7 ] * u.arcsec), 2 * u.arcsec,  0.2 * u.T ]

# To generate and view:
arr_Data = generate_example_data(arr_grid_shape, qua_xrange, qua_yrange, arrA2, arrA3)
aMap = dummyDataToMap(arr_Data, qua_xrange, qua_yrange)
aMap.peek()

auto_examplesimagessphx_glr_plot_gaussian_example_data_003.png

You can add as many poles as you want:

arr_Data = generate_example_data(arr_grid_shape, qua_xrange, qua_yrange, arrA0, arrA1, arrA2, arrA3)
aMap = dummyDataToMap(arr_Data, qua_xrange, qua_yrange)
aMap.peek()

auto_examplesimagessphx_glr_plot_gaussian_example_data_004.png

And being a map you can use all the normal SunPy functions, such as saving the map using aMap.save(filepath).

Total running time of the script: (0 minutes 0.911 seconds)

Download Python source code: plot_gaussian_example_data.py

Download IPython notebook: plot_gaussian_example_data.ipynb

Generated Data Extrapolation

Generated Data Extrapolation¶

In this example you will be generating some example data and extrapolate this using the basic potential extrapolator.

You can start by importing the necessary module components.

# Module imports
from solarbextrapolation.map3dclasses import Map3D
#from solarbextrapolation.potential_field_extrapolator import PotentialExtrapolator
from solarbextrapolation.extrapolators import PotentialExtrapolator
from solarbextrapolation.example_data_generator import generate_example_data, dummyDataToMap
from solarbextrapolation.visualisation_functions import visualise

You also need the ability to convert astropyunits and use MayaVi for visualisation.

# General imports
import astropy.units as u
from mayavi import mlab
import numpy as np

You are going to try and define a 3D cuboid grid of 20x22x20 with ranges in arcseconds, these parameters can be stored in the following lists and astropy quantities.

# Input parameters:
arr_grid_shape = [ 20, 22, 20 ]         # [ y-size, x-size ]
xrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
yrange = u.Quantity([ -11.0, 11.0 ] * u.arcsec)
zrange = u.Quantity([ 0,     20.0 ] * u.arcsec)

The generated data will consist of a 2D space with 2 Gaussian spots, one positive and one negative, on a background of 0.0. solarbextrapolation.example_data_generator provides many ways to achieve this, including letting it randomly generate the position, magnitude and size of each spot. In this case you will manually define the parameters of each spot as a list, using percentage units so that the spots will be inside the given ranges of any generated data:

# Manual Pole Details
#arrA# = [ position, size, maximum strength ]
arrA0 = [ u.Quantity([ 25, 25 ] * u.percent), 10.0 * u.percent,  0.2 * u.T ]
arrA1 = [ u.Quantity([ 75, 75 ] * u.percent), 10.0 * u.percent, -0.2 * u.T ]

You generate the data using generate_example_data(…) and create a map with this using dummyDataToMap(…).

# Generate the data and make into a map
arr_data = generate_example_data(arr_grid_shape[0:2], xrange, yrange, arrA0, arrA1)
map_boundary = dummyDataToMap(arr_data, xrange, yrange)

You can check the resulting generated data by using peek().

map_boundary.peek()

auto_examplesimagessphx_glr_plot_potential_extrapolation_of_example_data_001.png

You now simply want to extrapolate using this boundary data, this is achieved by first creating a potential extrapolator object and then by running the extrapolate on this to return a Map3D object with the resulting vector field.

# Use potential extrapolator to generate field
aPotExt = PotentialExtrapolator(map_boundary, zshape=arr_grid_shape[2], zrange=zrange)
aMap3D  = aPotExt.extrapolate(enable_numba=True)

# The Extrapolations run time is stored in the meta
floSeconds = np.round(aMap3D.meta['extrapolator_duration'],3)
print('\nextrapolation duration: ' + str(floSeconds) + ' s\n')

Out:

  True

extrapolation duration: 1.657 s

Note that you used enable_numba=True to speed up the computation on systems with Anaconda numba installed.

You can now get a quick and easy visualisation using the solarbextrapolation.example_data_generator.visualise tools:

# Visualise the 3D vector field
fig = visualise(aMap3D,
                boundary=map_boundary,
                volume_units=[1.0*u.arcsec, 1.0*u.arcsec, 1.0*u.Mm],
                show_boundary_axes=False,
                boundary_units=[1.0*u.arcsec, 1.0*u.arcsec],
                show_volume_axes=True,
                debug=False)
mlab.show()

auto_examplesimagessphx_glr_plot_potential_extrapolation_of_example_data_000.png

Out:

shape: (20L, 22L, 20L, 3L)

Note that the parameters here are simply to decide what boundary ranges to display.

Total running time of the script: (0 minutes 6.418 seconds)

Download Python source code: plot_potential_extrapolation_of_example_data.py

Download IPython notebook: plot_potential_extrapolation_of_example_data.ipynb

Potential Extrapolation from an Image

Potential Extrapolation from an Image¶

This example is to demonstrate using the potential extrapolator on an image. It was built for a bit of fun.

Traceback (most recent call last):
  File "/opt/miniconda/envs/solarb/lib/python2.7/site-packages/sphinx_gallery/gen_rst.py", line 467, in execute_script
    exec(code_block, example_globals)
  File "<string>", line 11, in <module>
ImportError: No module named potential_field_extrapolator

from __future__ import print_function

# General imports
from astropy import units as u
from scipy import misc
from mayavi import mlab
import numpy as np

# Module imports
from solarbextrapolation.potential_field_extrapolator import PotentialExtrapolator
from solarbextrapolation.example_data_generator import generate_example_data, dummyDataToMap
from solarbextrapolation.visualisation_functions import visualise

# The input parameters:
arr_grid_shape = [ 50, 50, 50 ]         # [ y-size, x-size ]
xrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
yrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
zrange = u.Quantity([ 0,     20.0 ] * u.arcsec)

# Manual Pole Details
arrA0 = [ u.Quantity([ 25, 25 ] * u.percent), 10.0 * u.percent,  0.2 * u.T ]
arrA1 = [ u.Quantity([ 75, 75 ] * u.percent), 10.0 * u.percent, -0.2 * u.T ]

# Generate the data and make into a map
arr_data = generate_example_data(arr_grid_shape[0:2], xrange, yrange, arrA0, arrA1)
print('\n' + str(type(arr_data.dtype)))
print(str(arr_data.shape))
print(str(arr_data.dtype) + '\n')
arr_image = (misc.imread('sunpy_powered_50x50.png')[...,:3] - 127.5) / 1270.0
arr_data = np.zeros(arr_image.shape[:2])
for i in range(0,arr_data.shape[0]):     # Row/Y
    for j in range(0,arr_data.shape[1]): # Column/X
        arr_data[i,j] = ((arr_image[i,j,0] + arr_image[i,j,1] + arr_image[i,j,2]) / 3.0)
print('\n' + str(type(arr_data.dtype)))
print(str(arr_data.shape))
print(str(arr_data.dtype) + '\n')

map_boundary = dummyDataToMap(arr_data, xrange, yrange)

# Use potential extrapolator to generate field
aPotExt = PotentialExtrapolator(map_boundary, zshape=arr_grid_shape[2], zrange=zrange)
aMap3D = aPotExt.extrapolate(enable_numba=True)
print('extrapolator_duration:' + str(aMap3D.meta['extrapolator_duration']))


# Visualise
visualise(aMap3D,
          boundary=map_boundary,
          volume_units=[1.0*u.arcsec, 1.0*u.arcsec, 1.0*u.Mm],
          show_boundary_axes=False,
          boundary_units=[1.0*u.arcsec, 1.0*u.arcsec],
          show_volume_axes=True,
          debug=False)
mlab.show()

Total running time of the script: (0 minutes 0.000 seconds)

Download Python source code: plot_potential_extrapolation_of_images.py

Download IPython notebook: plot_potential_extrapolation_of_images.ipynb

Defining and Run a Custom Analytical Model

Defining and Run a Custom Analytical Model¶

Here you will be creating trivial analytical model following the API.

You can start by importing the necessary module components.

# Module imports
from solarbextrapolation.map3dclasses import Map3D
from solarbextrapolation.analyticalmodels import AnalyticalModel
from solarbextrapolation.visualisation_functions import visualise

You also need the ability to convert astropyunits, manipulate numpy arrays and use MayaVi for visualisation.

# General imports
import astropy.units as u
import numpy as np
from mayavi import mlab

You are going to try and define a 3D cuboid grid of 20x22x20 with ranges in arcseconds, these parameters can be stored in the following lists and astropy quantities.

# Input parameters:
qua_shape   = u.Quantity([   20,   20,  20] * u.pixel)
qua_x_range = u.Quantity([ -80.0,  80 ] * u.Mm)
qua_y_range = u.Quantity([ -80.0,  80 ] * u.Mm)
qua_z_range = u.Quantity([   0.0, 120 ] * u.Mm)

From the above parameters you can derive the grid step size and total size in each dimension.

"""
# Derived parameters (make SI where applicable)
x_0 = x_range[0].to(u.m).value
Dx = (( x_range[1] - x_range[0] ) / ( tup_shape[0] * 1.0 )).to(u.m).value
x_size = Dx * tup_shape[0]
y_0 = y_range[0].to(u.m).value
Dy = (( y_range[1] - y_range[0] ) / ( tup_shape[1] * 1.0 )).to(u.m).value
y_size = Dy * tup_shape[1]
z_0 = z_range[0].to(u.m).value
Dz = (( z_range[1] - z_range[0] ) / ( tup_shape[2] * 1.0 )).to(u.m).value
z_size = Dy * tup_shape[2]
"""

You can define this analytical model as a child of the AnalyticalModel class.

class AnaOnes(AnalyticalModel):
    def __init__(self, **kwargs):
        super(AnaOnes, self).__init__(**kwargs)

    def _generate_field(self, **kwargs):
        # Adding in custom parameters to the metadata
        self.meta['analytical_model_routine'] = 'Ones Model'

        # Generate a trivial field and return (X,Y,Z,Vec)
        arr_4d = np.ones(self.shape.value.tolist() + [3])
        self.field = arr_4d

        # Extract the LoS Magnetogram from this:
        self.magnetogram.data = arr_4d[:,:,0,2]

        # Now return the vector field.
        return Map3D( arr_4d, self.meta )

You can instansiate a copy of the new analytical model.

aAnaMod = AnaOnes(shape=qua_shape, xrange=qua_x_range, yrange=qua_y_range, zrange=qua_z_range)

Note: you could use default ranges and grid shape using aAnaMod = AnaOnes().

You can now calculate the vector field.

aMap3D = aAnaMod.generate()

You can now see the 2D boundary data used for extrapolation.

aMap2D = aAnaMod.to_los_magnetogram()
aMap2D.peek()

auto_examplesimagessphx_glr_plot_define_and_run_trivial_analytical_model_001.png

You also visulise the 3D vector field:

fig = visualise(aMap3D,
                show_boundary_axes=False,
                show_volume_axes=False,
                debug=False)
mlab.show()

# Note: you can add boundary axes using:
"""
fig = visualise(aMap3D,
                show_boundary_axes=False,
                boundary_units=[1.0*u.arcsec, 1.0*u.arcsec],
                show_volume_axes=True,
                debug=False)
"""

auto_examplesimagessphx_glr_plot_define_and_run_trivial_analytical_model_000.png

Out:

shape: (20L, 20L, 20L, 3L)

Total running time of the script: (0 minutes 3.627 seconds)

Download Python source code: plot_define_and_run_trivial_analytical_model.py

Download IPython notebook: plot_define_and_run_trivial_analytical_model.ipynb

HMI FITS File Extrapolation

HMI FITS File Extrapolation¶

Example of extrapolating from a HMI fitts file using the potential extrapolator and visualising.

Traceback (most recent call last):
  File "C:\Users\alex_\Anaconda\lib\site-packages\sphinx_gallery\gen_rst.py", line 467, in execute_script
    exec(code_block, example_globals)
  File "<string>", line 25, in <module>
  File "C:\Users\alex_\Anaconda\lib\site-packages\sunpy-0.6-py2.7.egg\sunpy\map\map_factory.py", line 239, in __call__
    data_header_pairs, already_maps = self._parse_args(*args, **kwargs)
  File "C:\Users\alex_\Anaconda\lib\site-packages\sunpy-0.6-py2.7.egg\sunpy\map\map_factory.py", line 202, in _parse_args
    raise ValueError("File not found or invalid input")
ValueError: File not found or invalid input

# General imports
import numpy as np
import sunpy.map as mp
from astropy import units as u
from mayavi import mlab
import os

# Module imports
from solarbextrapolation.map3dclasses import Map3D
from solarbextrapolation.extrapolators import PotentialExtrapolator
from solarbextrapolation.visualisation_functions import visualise

# Cropping into the active region within the HMI map
str_vol_filepath = 'C:\\git\\solarbextrapolation\\examples\\2011-02-14__20-35-25__02_Bxyz.npy'
xrange = u.Quantity([50, 300] * u.arcsec)
yrange = u.Quantity([-350, -100] * u.arcsec)
zrange = u.Quantity([0, 250] * u.arcsec)
xrangeextended = u.Quantity([xrange.value[0] - 50, xrange.value[1] + 50] *
                            xrange.unit)
yrangeextended = u.Quantity([yrange.value[0] - 50, yrange.value[1] + 50] *
                            yrange.unit)

# Open the map and create a cropped version for the extrapolation.
map_hmi = mp.Map(
    'C:\\git\\solarbextrapolation\\examples\\2011-02-14__20-35-25__01_hmi.fits')
map_hmi_cropped = map_hmi.submap(xrange, yrange)
dimensions = u.Quantity([100, 100] * u.pixel)
map_hmi_cropped_resampled = map_hmi_cropped.resample(dimensions,
                                                     method='linear')

# Open the map and create a cropped version for the visualisation.
#map_boundary = mp.Map('C:\\git\\solarbextrapolation\\examples\\2011-02-14__20-35-25__02_aia.fits') # For AIA
map_boundary = mp.Map(
    'C:\\git\\solarbextrapolation\\examples\\2011-02-14__20-35-25__01_hmi.fits'
)  # For HMI

map_boundary_cropped = map_boundary.submap(xrangeextended, yrangeextended)

# Only extrapolate if we don't have a saved version
if not os.path.isfile(str_vol_filepath):
    aPotExt = PotentialExtrapolator(map_hmi_cropped_resampled,
                                    filepath=str_vol_filepath,
                                    zshape=dimensions[0].value,
                                    zrange=zrange)
    aMap3D = aPotExt.extrapolate()
aMap3D = Map3D.load(str_vol_filepath)
print('\nextrapolation duration: ' + str(np.round(aMap3D.meta['extrapolator_duration'], 3)) + ' s\n')

# Visualise this
visualise(aMap3D,
          boundary=map_boundary_cropped,
          scale=1.0 * u.Mm,
          boundary_unit=1.0 * u.arcsec,
          show_boundary_axes=False,
          show_volume_axes=True,
          debug=False)
mlab.show()

Total running time of the script: (0 minutes 0.000 seconds)

Download Python source code: plot_potential_extrapolation_of_hmi_fits_file.py

Download IPython notebook: plot_potential_extrapolation_of_hmi_fits_file.ipynb

Extrapolator Perfomance Testing

Extrapolator Perfomance Testing¶

In this example you will be running the potential field extrapolator both with numba enabled and disabled over a number of datasets and tabulating the results into an astropy table.

Note: if you don’t have conda numba installed the code should work but the results should not show any speed difference.

You can start by importing the necessary module components.

# Module imports
from solarbextrapolation.extrapolators import PotentialExtrapolator
from solarbextrapolation.example_data_generator import generate_example_data, dummyDataToMap

You also need the ability to convert astropyunits, use numpy arrays and astropy tables.

# General imports
from astropy import units as u
from astropy.table import Table
import numpy as np

You are going to create a series of volume grids with a given shape and then attribute arbitrary units to it’s axes.

lis_grid_shapes = [ [ 50, 50, 50 ] ]
xrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
yrange = u.Quantity([ -10.0, 10.0 ] * u.arcsec)
zrange = u.Quantity([ 0,     20.0 ] * u.arcsec)

Note that you could easily choose any grid dimensions:: e.g. [ [ 100, 100, 200 ] ]
or add more then one grid shape in the list:: e.g. [ [ 10, 10, 10 ],[ 50, 50, 50 ], [ 100, 100, 100 ] ]

to make the test more grid-size agnostic, but this will notably increase runtime.

You will create an example dataset using Gaussian spots, as show in the Generating Example Data example. We use the parameters:

# Manual Pole Details
arrA0 = [ u.Quantity([ 25, 25 ] * u.percent), 10.0 * u.percent,  0.2 * u.T ]
arrA1 = [ u.Quantity([ 75, 75 ] * u.percent), 10.0 * u.percent, -0.2 * u.T ]
arrA2 = [ u.Quantity([ 25, 75 ] * u.percent), 10.0 * u.percent,  0.1 * u.T ]
arrA3 = [ u.Quantity([ 75, 25 ] * u.percent), 10.0 * u.percent, -0.1 * u.T ]

# Generate the datasets and maps
#lis_maps = []
#lis_extrapolators = []

You will create an astropy table to store the runtimes of the extrapolations.

# A table for storing the data
t = Table(names=('grid size', 'time (min)', 'time (ave)', 'time (std)'), meta={'name': 'times tables'}, dtype=('S24', 'f8', 'f8', 'f8'))
t['time (min)'].unit = u.s
t['time (ave)'].unit = u.s
t['time (std)'].unit = u.s

You will store all the datasets to test with in a list. In this case the datasets will simply be the various generated example boundary data maps for the list of grid sizes, which is simply one example.

lis_datasets = []
for shape in lis_grid_shapes:
    lis_datasets.append([ str(shape), shape[2], zrange,
                          dummyDataToMap(generate_example_data(shape[0:2], xrange, yrange, arrA0, arrA1, arrA2, arrA3), xrange, yrange) ])

You may wish to run each test more than once, so you can use a parameter to autimate this.

int_trials = 1 # The times to repeat each extrapolation.

You iterate through the extrapolations on each dataset, adding teh runtime to the table.

for extrapolation in lis_datasets:
    # Setup the extrapolator and table
    aPotExt = PotentialExtrapolator(extrapolation[3], zshape=extrapolation[1], zrange=extrapolation[2])

    # List to store the trial
    lis_times = []

    # Run the extrapolation without numba for each dataset (map and ranges).
    for i in range(0, int_trials):
        aMap3D = aPotExt.extrapolate(enable_numba=False)
        lis_times.append(aMap3D.meta['extrapolator_duration'])
    t.add_row([extrapolation[0], np.round(np.min(lis_times), 2), np.round(np.average(lis_times), 2), np.round(np.std(lis_times), 2)])

    # List to store the trial
    lis_times = []

    # Run the extrapolation with numba for each dataset (map and ranges).
    for i in range(0, int_trials):
        aMap3D = aPotExt.extrapolate(enable_numba=True)
        lis_times.append(aMap3D.meta['extrapolator_duration'])
    t.add_row(['(numba)'+extrapolation[0], np.round(np.min(lis_times), 2), np.round(np.average(lis_times), 2), np.round(np.std(lis_times), 2)])

You can now see the results in the table.

print t

Total running time of the script: (0 minutes 0.000 seconds)

Download Python source code: potential_extrapolation_performance_tests.py

Download IPython notebook: potential_extrapolation_performance_tests.ipynb

Extrapolating HMI Data

Extrapolating HMI Data¶

In this example you will be downloading boundary data from VSO, extrapolating using the potential extrapolator and visualising in MayaVi.

You start by importing the necessary modules.

# General imports
import numpy as np
import sunpy.map as mp
from sunpy.net import vso
from astropy import units as u
from mayavi import mlab # Necessary for visulisation
import os

# Module imports
from solarbextrapolation.map3dclasses import Map3D
from solarbextrapolation.extrapolators import PotentialExtrapolator
from solarbextrapolation.visualisation_functions import visualise

You will retrieve the boundary data from the VSO using the SunPy VSO client. In this case we will retrieve an SDO HMI line-of-sight magnetogram that was made on the 14th of February 2011, as used in Sun et al (2012).

# Create a new VSOClient instance
client = vso.VSOClient()

# Build the query, this can return one item, or a list of them to DL matching
# the given filters.
result_hmi = client.query(
    # The following are filters for collecting the desired data.
    vso.attrs.Time((2011, 2, 14, 20, 34, 0), (2011, 2, 14, 21, 0, 0)), # Time range.
    vso.attrs.Instrument('HMI'),               # Helioseismic and Magnetic Imager.
    vso.attrs.Physobs('LOS_magnetic_field'),   # Physical observables
    vso.attrs.Sample(4000 * u.s)               # Only take a shot every $var seconds.
    # More observables at http://sdac.virtualsolar.org/cgi/show_details?keyword=PHYSOBS
)

# Save the results to fits files. (Using Rice compression if possible)
data_hmi = client.get(result_hmi, methods=('URL-FILE_Rice', 'URL-FILE')).wait()

You may also decide to get the corrisponding SDO AIA data showing the EUV image at the same time, this can be used to see the flux tubes for comparrison to the vector field streamlines for visulisation.

# Query VSO.
result_aia = client.query(
    vso.attrs.Time((2011, 2, 14, 20, 34, 0), (2011, 2, 14, 21, 0, 0)), # Time range.
    vso.attrs.Instrument('AIA'),    # Helioseismic and Magnetic Imager.
    vso.attrs.Physobs('intensity'), # Physical observables
    vso.attrs.Sample(4000 * u.s)    # Only take a shot every $var seconds.
    # More observables at http://sdac.virtualsolar.org/cgi/show_details?keyword=PHYSOBS
)

# Save the results to fits files. (Using Rice compression if possible)
data_aia = client.get(result_aia, methods=('URL-FILE_Rice', 'URL-FILE')).wait()

You want to crop on solar-x and solar-y the the active region of interest. Likewise you want to decide on the altertude ranges to extrapolate within. Extrapolators use astropy quantities for ranges, importanmtly these are designed to work with either physical length or angular units, conversion is done using the assumption the boundary data ios on the surface of the sun and following the small angle approximation. In this case we use angular uniits (arcsec specifically) for the zrange quantity, this is physically meaningless, but gives an easy way to ensure your zrange is similar to teh other ranges. We also want extended solar-x and solar-y ranges for plotting the

# Cropping into the active region within the HMI map
xrange = u.Quantity([50,    300] * u.arcsec)
yrange = u.Quantity([-350, -100] * u.arcsec)
zrange = u.Quantity([0,     250] * u.arcsec)

# Open the map and create a cropped version for the extrapolation.
map_hmi = mp.Map(data_hmi[0])
map_hmi_cropped = map_hmi.submap(xrange, yrange)

If your boundary data has a high resolution then you may need to resample to ensure it extrapolates within a reasonable timeframe.

# Resample boundary data map
shape = u.Quantity([20, 20] * u.pixel)
map_hmi_cropped_resampled = map_hmi_cropped.resample(shape, method='linear')

You can check the resulting generated data by using peek().

map_hmi_cropped_resampled.peek()

_images/sphx_glr_plot_potential_extrapolation_of_hmi_data_001.png

To speed up repeat usage of this script it will save the extrapolation output, you can use os.path.isfile() to check if the file already exists, assuming it doesn’t you will extrapolate and create it, otherwise you load it.

# Only extrapolate if we don't have a saved version
str_vol_filepath = data_hmi[0][0:-5] + '_Bxyz.npy'
if not os.path.isfile(str_vol_filepath):
    # Create the potential extrapolator and run the extrapolate method.
    aPotExt = PotentialExtrapolator(map_hmi_cropped_resampled, filepath=str_vol_filepath, zshape=20, zrange=zrange)
    aMap3D = aPotExt.extrapolate()
# Load the results.
aMap3D = Map3D.load(str_vol_filepath)
#print '\nextrapolation duration: ' + str(np.round(aMap3D.meta['extrapolator_duration'],3)) + ' s\n'

Out:

False

For the perposes of visualisation we will want an extended boundary data, not just that of the extrapolated region, and at the instruments full resolution, not resampled.

xrangeextended = u.Quantity([ xrange.value[0] - 50, xrange.value[1] + 50 ] * xrange.unit)
yrangeextended = u.Quantity([ yrange.value[0] - 50, yrange.value[1] + 50 ] * yrange.unit)

# Open the map and create a cropped version for the visualisation.
map_boundary = mp.Map(data_hmi[0])
map_boundary_cropped = map_boundary.submap(xrangeextended, yrangeextended)

You can now get a quick and easy visualisation using the solarbextrapolation.example_data_generator.visualise tools:

# Visualise the 3D vector field
visualise(aMap3D, boundary=map_boundary_cropped, scale=1.0*u.Mm, boundary_unit=1.0*u.arcsec, show_boundary_axes=False, show_volume_axes=True, debug=False)
mlab.show()

_images/sphx_glr_plot_potential_extrapolation_of_hmi_data_000.png

Total running time of the script: (0 minutes 54.427 seconds)

Download Python source code: plot_potential_extrapolation_of_hmi_data.py

Download IPython notebook: plot_potential_extrapolation_of_hmi_data.ipynb

Defining a Custom Preprocessor and Extrapolator

Defining a Custom Preprocessor and Extrapolator¶

Here you will be creating trivial preprocessor and and exztrqapolatoirs following the API.

You start by importing the necessary modules.

# General imports
import sunpy.map as mp
import numpy as np
from mayavi import mlab # Necessary for visulisation

# Module imports
from solarbextrapolation.preprocessors import Preprocessors
from solarbextrapolation.extrapolators import Extrapolators
from solarbextrapolation.map3dclasses import Map3D
from solarbextrapolation.visualisation_functions import visualise

Preprocessor Defining a trivial preprocessor that returns a zeros map for any given input map.

   class PreZeros(Preprocessors):
       def __init__(self, map_magnetogram):
           super(PreZeros, self).__init__(map_magnetogram)

       def _preprocessor(self):
           # Adding in custom parameters to the meta
           self.meta['preprocessor_routine'] = 'Zeros Preprocessor'

           # Creating the trivial zeros map of the same shape as the input map
           map_output = mp.Map((np.zeros(self.map_input.data.shape),
                                       self.meta))

           # Outputting the map.
           return map_output

Make an input map that we will run the preprocessor on.
This will be changed to using the sample HMI image.

aMap2D = mp.Map(‘C://git//solarextrapolation//solarextrapolation//data//example_data_(100x100)__01_hmi.fits’)

from solarbextrapolation.example_data_generator import generate_example_data, dummyDataToMap
import astropy.units as u
aMap2D = arr_Data = dummyDataToMap(generate_example_data([ 20, 20 ],u.Quantity([ -10.0, 10.0 ] * u.arcsec),u.Quantity([ -10.0, 10.0 ] * u.arcsec)), u.Quantity([ -10.0, 10.0 ] * u.arcsec), u.Quantity([ -10.0, 10.0 ] * u.arcsec))

Instansiate the preprocessor and process the input map.

aPrePro = PreZeros(aMap2D.submap([0, 10]*u.arcsec, [0, 10]*u.arcsec))
aPreProMap = aPrePro.preprocess()

You can plot the preprocessed map using peek.

aPreProMap.peek()

auto_examplesimagessphx_glr_plot_define_and_run_trivial_preprocessor_and_extrapolator_001.png

You can also access the metadata of the preprocessor like any map:

print "preprocessor_routine: " + str(aPreProMap.meta['preprocessor_routine'])
print "preprocessor_duration: " + str(aPreProMap.meta['preprocessor_duration'])

Traceback (most recent call last):
  File "C:\Users\alex_\Anaconda\lib\site-packages\sphinx_gallery\gen_rst.py", line 467, in execute_script
    exec(code_block, example_globals)
  File "<string>", line 1
    print "preprocessor_routine: " + str(aPreProMap.meta['preprocessor_routine'])
                                 ^
SyntaxError: invalid syntax

Extrapolator Defining a trivial extrapolator that returns a volume of one vectors.

class ExtOnes(Extrapolators):
    def __init__(self, map_magnetogram, **kwargs):
        super(ExtOnes, self).__init__(map_magnetogram, **kwargs)

    def _extrapolation(self):
        # Adding in custom parameters to the meta
        self.meta['extrapolator_routine'] = 'Ones Extrapolator'

        #arr_4d = np.ones([self.map_boundary_data.data.shape[0], self.map_boundary_data.data.shape[0], self.z, 3])
        arr_4d = np.ones(self.shape.tolist() + [3])
        return Map3D(arr_4d, self.meta)

Instansiate the preprocessor and extrapolate.

aExt = ExtOnes(aPreProMap, zshape=10)
aMap3D = aExt.extrapolate()

You can visulise the field using MayaVi.

fig = visualise(aMap3D,
                boundary=aPreProMap,
                show_boundary_axes=False,
                show_volume_axes=False,
                debug=False)
mlab.show()

"""

# aPreProData = aMap2D.submap([0,10], [0,10])

# Some checks:
#aPreProData.data # Should be a 2D zeros array.
#aPreProData.meta
#aPreProData.meta['preprocessor_routine']
#aPreProData.meta['preprocessor_start_time']

Traceback (most recent call last):
  File "C:\Users\alex_\Anaconda\lib\site-packages\sphinx_gallery\gen_rst.py", line 467, in execute_script
    exec(code_block, example_globals)
  File "<string>", line 17
    """

# aPreProData = aMap2D.submap([0,10], [0,10])

# Some checks:
#aPreProData.data # Should be a 2D zeros array.
#aPreProData.meta
#aPreProData.meta['preprocessor_routine']
#aPreProData.meta['preprocessor_start_time']

                                           ^
SyntaxError: EOF while scanning triple-quoted string literal

Testing an extrapolator

# Define trivial extrapolator
class ExtZeros(Extrapolators):
    def __init__(self, map_magnetogram, **kwargs):
        super(ExtZeros, self).__init__(map_magnetogram, **kwargs)

    def _extrapolation(self):
        # Adding in custom parameters to the meta
        self.meta['extrapolator_routine'] = 'Zeros Extrapolator'

        arr_4d = np.zeros([self.map_boundary_data.data.shape[0],
                           self.map_boundary_data.data.shape[0], self.z, 3])
        return Map3D((arr_4d, self.meta))


aExt = ExtZeros(
    aPreProData,
    filepath='C://Users/Alex/solarextrapolation/solarextrapolation/3Dmap.m3d')
aMap3D = aExt.extrapolate()

# Some checks:
#aMap3D.data # Should be a 4D zeros array.
#aMap3D.meta
#aMap3D.meta['extrapolator_routine']
#aMap3D.meta['extrapolator_start_time']

# Testing a Map3DCube

aMapCube = Map3DCube(aMap3D, aMap3D)
aMapCube[0]
aMapCube[0].data
aMapCube[0].meta
aMapCube[1].data
aMapCube[1].meta
"""

Traceback (most recent call last):
  File "C:\Users\alex_\Anaconda\lib\site-packages\sphinx_gallery\gen_rst.py", line 467, in execute_script
    exec(code_block, example_globals)
  File "<string>", line 36
    """
      ^
SyntaxError: EOF while scanning triple-quoted string literal

Total running time of the script: (0 minutes 0.261 seconds)

Download Python source code: plot_define_and_run_trivial_preprocessor_and_extrapolator.py

Download IPython notebook: plot_define_and_run_trivial_preprocessor_and_extrapolator.ipynb

solarbextrapolation.extrapolators Package¶

Classes¶

`Extrapolators`(map_magnetogram, **kwargs)	Common class for all 3D vector field extrapolation routines.
`PotentialExtrapolator`(map_magnetogram, **kwargs)	This is a greens function for extrapolating the potential (scalar) field above a given magnetogram.

Class Inheritance Diagram¶

solarbextrapolation.preprocessors Package¶

Classes¶

Preprocessors(map_data, **kwargs) A common class for all 2D pre-processing routines, tools used to pre-process the 2D sunpy map data for use in extrapolations.

Class Inheritance Diagram¶

solarbextrapolation.analyticalmodels Package¶

Classes¶

AnalyticalModel(**kwargs) Common class for the development of anylitical models of magnetic fields.

Class Inheritance Diagram¶

solarbextrapolation.utilities Module¶

Functions¶

`decompose_ang_len`(qua_input, **kwargs)	Function to help decompose quantities that have an equivilence between angles and length, such as photospheric observational angles and object sizes.
`si_this_map_OLD`(map)	Basic function to create a deep copy of a map but with all units in SI.
`si_this_map`(map)	Basic function to create a deep copy of a map but with all units in SI.

solarbextrapolation.map3dclasses Module¶

Classes¶

`Map3D`(data, meta, **kwargs)	A basic data structure for holding a 3D numpy array of floats or 3-float vectors and metadata.
`Map3DCube`(args, *kwargs)	A basic data structure for holding a list of Map3D objects.
`Map3DComparer`(map3D, args, *kwargs)	Class for comparrison of vector fields.