tsssss / geopack Goto Github PK

Hi I'm having some problems using some of the functions in geopack, for example

import geopack as gp

xgsm, ygsm, zgsm = gp.geogsm(x,y,100,1)
print(xgsm,ygsm,zgsm)

returns module 'geopack' has no attribute 'geogsm' as well as other functions such as igrf_geo similarly saying that geopack doesn't have an attribute
I used python3 -m pip install geopack to install it, is there anything else I need to install for geopack to work.

thanks in advance.

Hi,

Geopack is incompatible with numpy after v1.24.0, since np.float or similar aliases are deprecated. Could you please make an update on this issue?

I typically refer to https://heliopython.org/projects/ when looking for space physics related python code. I found out about the existence of this code by reading someone else's code. I think this is a very useful project and the community would benefit if more people knew about it.

The IDL and Python outputs from the T01 model are differing by a few nanotesla for a large part of the 2007-03-23 THEMIS-A ephemeris which I'm testing with. The differences are occurring in regions with low field strength (sometimes 16x the reference field value), so I think this is more than just roundoff error.

See the attached plot: top two panels are the IDL and Python T01 field outputs; third panel is the difference, bottom panel shows the GSM positions in RE used as model inputs. The other parameters were: pdyn = 2.0, dsti = -30.0, yimf = 0.0, zimf = -5.0, g1 = 6.0, g2 = 10.0.

Hello!

I ran into some differences with T96 while comparing the results to IDL; to fix the problems, I updated to the GitHub release with:

pip install https://github.com/tsssss/geopack/archive/master.zip

instead of using:

pip install geopack

So I think the last T96 bug fix hasn't made it to a release on PyPI yet.

Are there plans to do a new release soon?

Thanks in advance!

Would it be possible to make geopack available on conda-forge? We would like to be able to use pyspedas on heliocloud, and it's a requirement that our package and all dependencies be available in conda.

I hear it's possible to add packages to conda-forge, even if you're not the package owner -- I am willing to try adding it myself, but only if that's OK with you.

Hello,
I am testing and comparing the outputs of this geopack for two different models with those of IDL and I am getting very different results.

I have attached the figures which I plotted for this purpose. For each panel in the figure I have plotted the difference between the output from IDL and Python (b_out_idl - b_out_py). Difference from T96 are plotted in the first row whereas the difference between T01 model are plotted in the second row. The value of x_gsm is at the top left position whereas the value of DST used for the code is at top right

For both T96 and T01 model far from the earth, into the tail, the difference between two outputs is minimal, almost always less than 2nT, and thus can be ignored in most cases (Figure 1, for x_gsm = 9.98 R_E).

However, there seem to be significant difference in output when we move close to the Earth or when we are on the dayside. Though the difference is smaller for T01 outputs compared to those of T96, they are still significant, specially far from the Earth on the dayside (Figure 2).

I wonder if anyone else have had similar issues.

The parameters I used for the codes are as follows:

time = 2015-01-01 00:00:00 GMT
par = [5, 0, 1, 1, 0, 0, 0, 0, 0, 0]
x_gsm = np.linspace(-15.1, 15, 61)
y_gsm = np.linspace(-15.1, 15, 61)
z_gsm = np.linspace(-15.1, 15, 61)

I can provide the full python and IDL codes if that will help with the reproduction of these differences.

Figure 1

Figure 2

Figure 3

Do you have any preferences on how this implementation of geopack should be referenced?

Please update to IGRF-13.

It would be great to have the ability to evaluate the Tsyganenko models at multiple points using vectorized processing. For example, evaluating B for every point on a grid. Some scientific use cases:

1. Model comparison in the global perspective
2. Evaluating the models on a grid, and then plugging the result into other code that uses gridded output

This should be possible if operations are all done in a way that supports arrays in addition to individual points. Right now, I can iterate on every point in the grid, but it's very slow to do it this way because the Python interpreter is doing a lot of unnecessary work.

Maybe there are some errors in t96.py?
one error with BIRK1TOT_02 fuction in t96.py, should change c2[54]=29.2043 to c2[54]=-29.2043;
one error with condip1(xi) fuction in t96.py,should change if (i == 2) | (i == 4) | (i == 6): to if (i == 2) | (i == 4) | (i == 5):

Hi,

I am trying to do a field line trace and get the points for the entire field line. However, I am only getting a single set of coordinates even though the documentation for the function says it should also output a vector of coordinates. Could you please help?

Example Code:
geopack.trace(1,0,.5, 1,10,1,2,'t89','igrf')

Output:
(0.9921003911113675, 0.00048052726154043516, 0.1254269006978363)

Thank you!

My install of geopack failed on a M1 mac because the scipy install failed.

I was able to get it to work using

pip install --pre -i https://pypi.anaconda.org/scipy-wheels-nightly/simple scipy

This isn't really a bug report but rather a comment to help others that may run into this issue.

In my code that uses geopack, I've added

import platform
import magnetovis as mvs

platform_str = platform.platform()
if platform_str.endswith('arm64-arm-64bit'):
    try:
        from geopack import geopack
    except:
        msg = 'pip install --pre -i https://pypi.anaconda.org/scipy-wheels-nightly/simple scipy'
        msg = '\n\nThis function uses a library that depends on scipy. On Mac M1, you need to install scipy using: \n\n   ' + msg + '\n'
        raise ImportError(msg)

One thing to consider is the only place scipy is used is in t96.py. There, only scipy.special is used for Bessel functions. The scipy dependency adds 40 MB to the required packages and so one may want to consider using https://libraries.io/pypi/special-functions, which is only 2 MB and has no external dependencies.

Hi, I am trying to use the geopack.trace function to calculate the limiting field line position antiparallel to my point (dir=1). However whenever I use dir=1 I get the following error:

UnboundLocalError                         Traceback (most recent call last)
<ipython-input-121-b09879ee0d08> in <module>
      1 parmod = [data['Pdyn'][k],data['SYMH'][k],data['BYGSM'][k],data['BZGSM'][k],data['W1'][k],data['W2'][k],data['W3'][k], \
      2           data['W4'][k],data['W5'][k],data['W6'][k]]
----> 3 a,b,c,_,_,_ = geopack.trace(xgsm,ygsm,zgsm,dir=1)
      4 a1,b1,c1 = geopack.geogsm(a,b,c,-1)
      5 print(a1,b1,c1)

~/.conda/envs/rtenv36/lib/python3.6/site-packages/geopack/geopack.py in trace(xi, yi, zi, dir, rlim, r0, parmod, exname, inname, maxloop)
   1165             # find the footpoint position by interpolating between the current and previous field line points:
   1166             r1=(r0-r)/(rr-r)
-> 1167             x=x-(x-xr)*r1
   1168             y=y-(y-yr)*r1
   1169             z=z-(z-zr)*r1

UnboundLocalError: local variable 'xr' referenced before assignment

When I look at the source code, I see that xr, yr, zr are all never defined but are referenced before their assignment. Has part of the source code been omitted?

Hi

I'm testing this python geopack package with the IDL routine.

I first tested with the parameters used in the test_geopack1.md as following:
Time/dipole tilt: 2001-01-01/02:03:04 UT/ps = -0.533585131 rad
Position: (xgsm,ygsm,zgsm) = (-5.1,0.3,2.8)
Model parameter:
T89: iopt = 2
T96, T01, and T04: par = [2., -87., 2., -5., 0., 0., ps, xgsm,ygsm,zgsm].

t89 | idl | 20.772 | -0.647 | -15.072
  | python | 20.772 | -0.647 | -15.072
t96 | idl | 61.179 | -1.461 | -40.451
  | python | 61.179 | -1.461 | -40.451
t01 | idl | 46.973 | 1.544 | -31.355
  | python | 46.354 | 1.440 | -32.000
t04 | idl | -8.080 | -1.512 | -0.722
  | python | 12.010 | -2.635 | -11.964

This result is similar to what is in test_geopack1.md except t04 for IDL routine.

I then tested for one of the MMS data, and the results are much different
datetime | 2016-01-01/12:13:20
iopt / kp | 3 / 2
(xgsm,ygsm,zgsm) = | 6.96, 0.49, -0.66
parmod | 0.976, -50.20, 6.59, -5.512, 0., 0., tilt, x,y,z

t89 | idl | -3.657066785 | 1.085712758 | 13.86057416
  | python | -2.37842015 | 1.038191716 | 14.38198039

t96 | idl | 22.97810632 | -1.106448782 | -4.147864164
  | python | 50.00624944 | 6.397909797 | -147.0215389

t01 | idl | 5.001211425 | 0.02607965332 | 3.264103442
  | python | 6.230733691 | 4.283414418 | -2.58591528

t04 | idl | -15.25291193 | -8.650237035 | -10.97633666
  | python | 7.905332371 | -1.161844446 | -7.933264503

Could you help me with it?

Thanks

Hello! I found that the func t96 of geopack python version has very different results with the IDL version at some places, such as dayside magnetosphere. But at other places, they get exactly the same results. I began to doubt some parts of the codes have some problems.

After long time of debugging, I found two bugs:

1. t96.py: func r2_birk: line 1410: "r2outer" should be "r2sheet".
2. t96.py: func birk1tot_02: line 844: "[x1,y1,z1,ps]" should be "[x2,y2,z2,ps]".

The IDL version of geopack is directly compiled from the original FORTRAN77 codes, so I think its results are correct.
I diecrtly compared the FORTRAN77 codes with the python codes, and found the two bugs above.

Hi,

During packaging https://github.com/spedas/pyspedas for Guix I've faced with the
issue of failing unit test of geopack which comes as install dependency of the
project:

Using pytest
============================= test session starts ==============================
platform linux -- Python 3.10.7, pytest-7.1.3, pluggy-1.0.0 -- /gnu/store/jh59fh48mcffyz5wmsjj0p96xkkflbz0-python-wrapper-3.10.7/bin/python
cachedir: .pytest_cache
hypothesis profile 'default' -> database=DirectoryBasedExampleDatabase('/tmp/guix-build-python-geopack-1.0.10.drv-0/geopack-1.0.10/.hypothesis/examples')
rootdir: /tmp/guix-build-python-geopack-1.0.10.drv-0/geopack-1.0.10
plugins: hypothesis-6.54.5
collecting ... collected 1 item

geopack/test_geopack1.py::KnownValues::test_to_known_values FAILED       [100%]

=================================== FAILURES ===================================
_______________________ KnownValues.test_to_known_values _______________________

self = <geopack.test_geopack1.KnownValues testMethod=test_to_known_values>

    def test_to_known_values(self):
        """geopack should give known result with known input"""

        # test recalc, which returns dipole tilt angle in rad.
>       self.assertEqual(self.test['recalc'], geopack.recalc(self.test['ut'], *self.test['v_gse']))
E       AssertionError: -0.53356312486214 != -0.5335631237417825

geopack/test_geopack1.py:62: AssertionError
=========================== short test summary info ============================
FAILED geopack/test_geopack1.py::KnownValues::test_to_known_values - Assertio...
============================== 1 failed in 0.75s ===============================

• Python 3.10.7

• [email protected]

• [email protected]

• [email protected]

• [email protected]

Hello! Could you do a release to get the fix in #14 into PyPI?