dstansby / pfsspy Goto Github PK

As spotted over in #211, if the input has any NaN values these get propagated all the way through the solution. This isn't very helpful, so instead the input should be checked and a ValueError raised if any of the input data is non-finite.

Add a function called load_map(source, dtime), that automatically downloads (via. SunPy), loads the data, and returns an Input pfsspy object. Sources should be:

• - HMI
• - GONG
• - ADAPT (HMI/GONG)
• - MDI
• - LMASL
• - WSO
• - MWO

Hello,
I found 2 errors in the example code namely "open_closed_map.ipynb".

1. ----> 2 tracer = tracing.FortranTracer()
   AttributeError: module 'pfsspy.tracing' has no attribute 'FortranTracer'.

2. ----> 9 pols = field_lines.polarities.reshape(2 * nsteps + 1, nsteps + 1).T
   NameError: name 'field_lines' is not defined

How to fix these errors please?

ADAPT maps come with a few different realisations of the model that produces them. I think these would be best stored in a MapSequence, but I don't know if it's possible to register a map type to the SunPy map factory that will then return a MapSequence.

Todo

• Try loading a ADAPT file with sunpy.map.Map. Document what goes wrong
• Work out if it's possible to create a sunpy map source that returns a MapSequence
• I suspect the answer to this is 'no', and that some very non-trival changes would need to be made to sunpy to support this. If this is the case, writing an example for the pfsspy example gallery showing how to use sunpy.io to read in an adapt fits file, and put it in a MapSequence would be a great outcome.

In theory MDI synoptic magnetograms might exist somewhere, but they don't seem to be on JSOC: http://jsoc.stanford.edu/MDI/MDI_Magnetograms.html. If anyone knows where to get them, please post below!

At the moment tracing is done on a cartesian grid that doesn't quite cover the whole spherical domain. An extra layer should be added to the cartesian grid to make it cover the whole domain.

I started a HMI searching/downloading example in #161. This needs to be expanded to include a brief discussion of the different synoptic maps available from HMI, and the different keys needed to search for them. Information available here: http://jsoc.stanford.edu/HMI/LOS_Synoptic_charts.html

If anyone else wants to do this, feel free to use my PR linked above as the base for a new PR.

Hello,
I am running the example codes of PFSS but I am getting this error:
`ModuleNotFoundError Traceback (most recent call last)
in
10 from pfsspy import coords
11 from pfsspy import tracing
---> 12 from gong_helpers import get_gong_map

ModuleNotFoundError: No module named 'gong_helpers'`
I have uninstalled PFSS and reinstalled it; however, the error is still there.
Can you please check it?

Currently calculating magnetic field streamlines is quite slow; it would be good to optimize the calculation, possibly including using cython?

It seems map.meta is modified in place when the map, or some of it's inbuilt methods are called, and the 'KEYCOMMENTS' attribute is dropped. This is inconvenient because GongSynopticMap uses this attribute in it's __init__() method so if we do :

>>> import sunpy.map
>>> import pfsspy
>>> gongmap = sunpy.map.Map("<path/to/gong.fits>")
>>>  gongmap.meta.get('KEYCOMMENTS','Does not exist') 
Out[15]: 
{'SIMPLE': 'Written by IDL:  Sun Oct 27 09:04:39 2019',
 'BITPIX': 'Bits per pixel',
 'NAXIS': 'Number of axes',
...
>>> ## So far so good, now try using a gongmap method :
>>> print(gongmap.bottom_left_coord)                                                             
<SkyCoord (HeliographicCarrington: obstime=2019-03-10T00:00:00.000): (lon, lat, radius) in (deg, deg, km)
    (333.5, -83.95715372, 695700.)>
>>> ## Try accessing 'KEYCOMMENTS' again :
>>> gongmap.meta.get('KEYCOMMENTS','Does not exist') 
Out[19]: 'Does not exist'
>>> # Now we are missing information we need to re-instantiate the class!
>>> gongmap_remake = sunpy.map.Map(gongmap.data,gongmap.header)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
...
KeyError: 'keycomments'

This makes modifying a GongSynopticMap e.g. as in #169 tricky.

I haven't figured out yet how or why sunpy.map throws away this info, but perhaps the path of least resistance would be to make the GongSynopticMap subclass independent of the 'keycomments' field?

In the example showing how to overplot field lines on an AIA image, the Tx and Ty values are pulled out and the plot command is used. Instead, the coordinates can be fed directly to the plot_coord command. See the sunpy-1.0-paper figure (which was adapted from this example).

It would be good to test against sunpy (and probably astropy) dev to catch any changes before they're released.

These should have differently named methods to get their footpoints.

Description

The readthedocs front page is well designed with all relevant information. The GitHub front page as rendered from the README.md file is very minimal, does not provide information for the installation and/or documentation link. I suggest having more information upfront on the GitHub page helps the users navigate to the documentation and installation.

This issue is a non-essential suggestion as part of my JOSS review.

I cloned pfsspy from the github landing page (master branch) and found the dtime keyword appears to be missing from the Input class init() function . Seems an older version of it may have propagated into the master branch? No problems in the release versions.

Could probably do with a dedicated installation page, that includes

• Only python >=3.5
• Only on pip

The outputted magnetic field values should have coordinate system information attached. This probably involves using astropy coordinates in a clever way to represent a 3D vector field.

Where we know the synoptic map was made from data observed from earth, we should set their observer coordinates.

This involves:

• Checking if streamline calculation can be done in parallel (I think that the interpolator might be non-re-entrant at the moment?)
• Accepting a n x 3 input to the tracing
• Splitting the input coordinates into sub-arrays of coordinates
• Handing the sub-arrays to multiprocessing
• Allowing specification of the number of processes

I'm finding in the current release that all field lines are either positive or negative polarity and no distinction is made between open and closed field lines. The closed field lines appear to have an arbitrary polarity but I'm not sure where it is generated from - there's probably an easy fix which was implemented in previous releases to check if the field line has two photospheric connections or not. The open field lines appear to have the correct polarities still. See the final images in the dipole source solution or the GONG PFSS extrapolation examples.

I'm finding for applications where I am running pfsspy in a loop, e.g. to iterate over many days of magnetograms, my memory is filling up steadily and eventually crashes.

I've found at least one repeatable test case where I generate a pfsspy.Output object and then access it's contents. The attached image shows some code I run on the left and my system monitor on the right. I first run a loop where I create the Output object, po, and assign the bc property to a new variable and repeat this 4 times, always overwriting the same variables. In the system monitor you can see 4 pulse of CPU activity, the first time the memory usage increases, the remaining 3 it stays more or less constant, which is expected.

However, I next do the same loop but instead of accessing bc, I access bg. This is the next 4 pulses of cpu activity and you can see the memory usage climbs each time even though the variables are being overwritten in the same way.

Looking at the pfsspy.Output source code, the difference seems to be that bg has a functools.lru_cache decorator. I don't fully understand what this does but it seems to be causing this memory "leakage". Is there a solution to prevent this? What is the function of this decorator?

Currently the SkyCoord input to tracing has to be a 1D coordinate array. It should be pretty easy to accept an arbitrary shaped array - flatten it on input, and then reshape the output to match the input shape.

Hello,
I found the following 2 errors in the example code namely "tracer_performance.ipynb":

1.   2 tracers = [pfsspy.tracing.PythonTracer(),
   

----> 3 pfsspy.tracing.FortranTracer()]
AttributeError: module 'pfsspy.tracing' has no attribute 'FortranTracer'

2. ----> 2 ax.scatter(nseeds[1:len(times[0])], times[0][1:], label='python')
   NameError: name 'nseeds' is not defined

How to fix these errors please?

I thought I would open a new issue to track the fixes we need to make to HMI maps. See #158 for an example of how to implement a map source.

Running the example code in #161 the errors I get are:

(dev) Davids-MBP-3:doc dstansby$ python ../examples/finding_data/plot_hmi.py 
WARNING: FITSFixedWarning: CRDER1 = 'nan ' 
a floating-point value was expected. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: CRDER2 = 'nan ' 
a floating-point value was expected. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'Degree' -> 'deg'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'celfix' made the change 'In CUNIT2 : Invalid symbol in EXPON context in 'Sine Latitude''. [astropy.wcs.wcs]
Traceback (most recent call last):
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 566, in _do_parse
    return cls._parse_unit(s, detailed_exception=False)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 487, in _parse_unit
    raise ValueError()
ValueError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 459, in _get_unit
    return cls._parse_unit(t.value)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 483, in _parse_unit
    raise ValueError(
ValueError: Degree is not a valid unit. Did you mean degree or exadegree?

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/core.py", line 1854, in __call__
    return f.parse(s)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 551, in parse
    result = cls._do_parse(s, debug=debug)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 569, in _do_parse
    return cls._parser.parse(s, lexer=cls._lexer, debug=debug)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/extern/ply/yacc.py", line 331, in parse
    return self.parseopt_notrack(input, lexer, debug, tracking, tokenfunc)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/extern/ply/yacc.py", line 1061, in parseopt_notrack
    lookahead = get_token()     # Get the next token
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/extern/ply/lex.py", line 350, in token
    newtok = func(tok)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 165, in t_UNIT
    t.value = cls._get_unit(t)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/format/generic.py", line 461, in _get_unit
    raise ValueError(
ValueError: At col 0, Degree is not a valid unit. Did you mean degree or exadegree?

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "../examples/finding_data/plot_hmi.py", line 24, in <module>
    hmi_map.plot()
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/decorators.py", line 234, in wrapper
    return_ = wrapped_function(*func_args, **func_kwargs)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/sunpy/map/mapbase.py", line 1719, in plot
    axes = wcsaxes_compat.gca_wcs(self.wcs)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/sunpy/map/mapbase.py", line 294, in wcs
    w2.wcs.cdelt = u.Quantity(self.scale)
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/sunpy/map/mapbase.py", line 706, in scale
    return SpatialPair(self.meta.get('cdelt1', 1.) * self.spatial_units[0] / u.pixel,
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/sunpy/map/mapbase.py", line 714, in spatial_units
    return SpatialPair(u.Unit(self.meta.get('cunit1')),
  File "/Users/dstansby/miniconda3/envs/dev/lib/python3.8/site-packages/astropy/units/core.py", line 1874, in __call__
    raise ValueError(msg)
ValueError: 'Degree' did not parse as unit: At col 0, Degree is not a valid unit. Did you mean degree or exadegree? If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see http://docs.astropy.org/en/latest/units/combining_and_defining.html

Currently field lines traced through the PFSS solution reach slightly outside the boundary, since they use finite difference steps. These should be clipped so they start at exactly 1 solar radius, and finish at exactly the source surface. This involves interpolating a point between the last two points of the field line.

See #164 for the relevant JSOC keyword.

Similar to #201 , I think it would be handy to be able to arbitrarily sample a point inside the model and get the vector B field. This is already doable with the built-in interpolator _brgi but I had to do some digging to understand the needed input shape and to understand the output is cartesian while the input is a spherical coordinate.

Perhaps there could be a new method for Output objects which takes an arbitrary SkyCoord position, checks it is inside the model and returns Br,Bth,Bph or Bx,By,Bz depending on a keyword argument?

Happy to pull a PR together if you think it would be a good addition @dstansby

I recently realized that you can add parallel=True to the numba.jit decorator to automatically take advantage of multiple threads in JIT-ed computations. If you have loops that use range/arange, you can then replace these with numba.prange to parallelize them. I recently applied this in my much more rudimentary extrapolator and got a factor of 4 speedup for a 100-by-100-by-100 extrapolation (with 36 threads). As pfsspy already uses numba on some functions with a few similar-looking loops, it may be worth experimenting to see how much of a difference this makes performance-wise.

I am not sure if this is already an open issue or not -- if so, I'll close the issue.

I can't seem to create a PFSS model, even after implementing the fixes in the sunpy/sunpy#4053. The model seems to run, but the output is empty. Here is my code. Apologies if I'm making an obvious error!

https://codecov.io/gh/dstansby/pfsspy/src/2cfc38e6884e0171e3a0a7fb88d6ea3e8bd17334/pfsspy/__init__.py shows that the numba lines aren't being tested. One of the tests should install numba.

When downloading HMI data, two files are downloaded:

['/home/circleci/sunpy/data/hmi.mrsynop_small_720s.2210.synopMr.fits', '/home/circleci/sunpy/data/hmi.mrsynop_small_720s.2210.epts.fits']

I'm not sure what the second one is. This issue is to find out, and add a quick description to the search for HMI data example.

Not quite sure what form these should take, but there should be some tests for the field tracing code.

Currently the test in pfsspy/tests/test_pfsspy.py does absolutely nothing. A simple smoke test should be added, probably based of the single example in the examples directory.

This will probably be a section on the front page that includes a small bit of citation text, and a link to a zenodo DOI (which needs creating).

Hello,
I would suggest adding functionality to get the 3 components of the extrapolated magnetic field along each generated field line i.e., B(x,y,z)
I think that would be very helpful to trace the evolution of the magnetic field at different locations and at different heights.

Thank you very much for this amazing work!

I think we're going to need a custom client to be written for GONG synoptic maps. @abhijeetmanhas I think you were thinking of working on this since you've done similar work in SunPy for non-synoptic GONG products (sunpy/sunpy#3811)?

This work should probably be done in https://github.com/sunpy/sunpy instead of pfsspy.

Currently the field line tracing doesn't go quite to the boundary. It's possible to iterate to smaller step sizes near the boundary to reach the boundary.

Hi David,

I really want to use the pfsspy package to overplot field lines on AIA images i.e. this example: https://pfsspy.readthedocs.io/en/stable/auto_examples/using_pfsspy/plot_aia_overplotting.html#sphx-glr-auto-examples-using-pfsspy-plot-aia-overplotting-py

I am having problems importing the pfsspy.tracing module. I get the following errors when trying to import:

`import pfsspy
import pfsspy.coords as coords
import pfsspy.tracing as tracing
from gong_helpers import get_gong_map
Traceback (most recent call last):

File "", line 3, in
import pfsspy.tracing as tracing

ModuleNotFoundError: No module named 'pfsspy.tracing'`

I also get a similar error message for the gong_helpers module.

I installed pfsspy using pip. I have gfortran and streamtracer installed and updated various python packages. Is this something you have come across before and any ideas on how to fix this? It might be something simple that I have missed...

Steph

Roll Function
Take sunpy.map.Map instance, verify is synoptic map, edit Map.data (inplace?) : roll carrington longitude = 0 to LH edge of map

Interpolate to sinlat function
Take sunpy.map.Map instance, verify is synoptic map, verify latitudinal binning is not sinlat, edit Map.data (inplace?) : interpolate data to uniform binning in sinlat (same dimensions as previous?)

It would be convenient if pfsspy.Input accepted a HGC SunPy map as the input rather than a NumPy array.

https://stackoverflow.com/questions/41220617/python-3d-interpolation-speedup has example code to use cython for interpolation, which should massively speed up field line tracing.

Hi Dave,

Good stuff! Not really being a python user yet, I just tried out the examples, but got an error message for each of the two in the last step :-(

I'm quite sure it's somehow my fault - can you advise? Thanks!

Daniel

--

Currently the readthedocs build is failing due to a lot of memory being taken up, I'm guessing by streamline tracing. It should be possible to do a memory profile of the tracing to try and reduce memory usage.