T89,01 model

[kmrs-iitd]

Dear Pyspedas Experts,
I'm attempting to run the T89 model using my data, but I’m encountering extremely high values for the calculated magnetic field (around 10^45 nT). I have checked my data and confirmed that there are no outliers, NaN, or null values. Can someone assist me please with this issue?

code:

import pyspedas
import numpy as np
from pytplot import store_data, tplot, join_vec
from pyspedas.geopack.get_tsy_params import get_tsy_params
from pyspedas import tinterpol
from pyspedas.geopack import tt01,tt89,tt96

https://drive.google.com/file/d/1TXL7Cu--ear8peVNb0yGPE0oeJGugReL/view?usp=sharing :: this data file is shared through google drive
T89=pd.read_csv('~/Documents/T89.csv')

data = np.column_stack((T89['XGSM'], T89['YGSM'], T89['ZGSM']))
store_data('X_Y_Z_gsm', data={'x': T89['Date'], 'y': data})

tt89('X_Y_Z_gsm')

tplot('X_Y_Z_gsm_bt89')

The T01 model is also encountering an issue and displaying the following error:

The name X_Y_Z_gsm-itrp_bt01 is currently not in pytplot
The name X_Y_Z_gsm-itrp_bt01 is currently not in pytplot
Variable not found: X_Y_Z_gsm-itrp_bt01

Code:
data = np.column_stack((T89['XGSM'], T89['YGSM'], T89['ZGSM']))
store_data('X_Y_Z_gsm', data={'x': T89['Date'], 'y': data})

join_vec(['BX_GSE','BY_GSM', 'BZ_GSM'], newname='BX_BY_BZ_GSM_joined')

params = get_tsy_params('kyoto_dst',
'BX_BY_BZ_GSM_joined',
'proton_density',
'flow_speed',
't01', # or 't96', 'ts04'
pressure_tvar='Pressure',
speed=True)

tinterpol('X_Y_Z_gsm', 'proton_density')

tt01('X_Y_Z_gsm-itrp', parmod=params)

tplot('X_Y_Z_gsm-itrp_bt01')

[jameswilburlewis]

When I've seen results like this in pyspedas, it's usually because I've specified the input positions in units of Re rather than km. This is inconsistent with the convention in IDL SPEDAS, and we're working to fix this in pyspedas. Try passing the positions in units of km and see if that works better. You can see detailed examples of how pyspedas geopack expects its arguments to be processed in the unit tests under pyspedas/geopack/tests. Hope that helps... let us know if you're still having problems after trying different units!

[kmrs-iitd]

Thank you! I’ll give your suggestions a try and will provide an update if the issue continues.

[kmrs-iitd]

Sorry for delay. Actually, I am attempting to compare the magnetic field variations at Earth's polar regions using different models, specifically T89, T96, T01, and TS04. I converted coordinates from [X, Y, Z]_gsm [Re] to [X, Y, Z]_gsm [km].
The T89 model is working well, but I am still encountering issues with the T01 model.
However, both the T96 and TS04 models are working but producing unusually high values for the magnetic field, as shown in the images.

The T01 model gives blank screen with:

tinterpol (linear) was applied to: Pressure_interp
G variables required for T01 model; create a tplot variable containing the G variables, and provide the name of that keyword to the g_variables keyword.
tinterpol (linear) was applied to: X_Y_Z_gsm-itrp
parmod keyword required.
That name is currently not in pytplot
That name is currently not in pytplot
Variable not found: X_Y_Z_gsm-itrp_bt01

I applied the same data to all of these models:

https://drive.google.com/file/d/1TXL7Cu--ear8peVNb0yGPE0oeJGugReL/view?usp=sharing
This data file is shared through google drive

df = pd.read_csv('~/Documents/T89.csv')

Code:

--Changing Re--> Km

df['XGSM1'] = df['XGSM']*6371.2
df['YGSM1'] = df['YGSM']*6371.2
df['ZGSM1'] = df['ZGSM']*6371.2

data = np.column_stack((df['XGSM1'], df['YGSM'1], df['ZGSM1']))

store_data('X_Y_Z_gsm', data = {'x': df['Date'], 'y': data})

join_vec(['BX_GSE','BY_GSM', 'BZ_GSM'], newname='BX_BY_BZ_GSM_joined')

params = get_tsy_params('kyoto_dst',
'BX_BY_BZ_GSM_joined',
'proton_density',
'flow_speed',
't01', # or 't96', 'ts04'
pressure_tvar='Pressure',
speed=True)

tinterpol('X_Y_Z_gsm', 'proton_density')

tt01('X_Y_Z_gsm-itrp', parmod=params)

tplot('X_Y_Z_gsm-itrp_bt01')

I used same code for T96 and TS04.

Please help.

[jameswilburlewis]

That's a lot of data to search through.....could you point me to a few specific points where the model outputs are not what you expect? That would make it easier for me to reproduce the issue, and see if there's a problem with any of the models.

[kmrs-iitd]

Apologies for the confusion. Currently, I am struggling with T01 model, that is not producing results. I've now filtered the data (from one year) to cover just one month.

df = pd.read_csv('~/Documents/T89.csv')

df1 = df[(df['Date'] >= '1997-09-01 00:00:00') & (df['Date'] <= '1997-10-01 00:00:00')]

df1 = df1.drop_duplicates(subset=['Date']).sort_values('Date').reset_index(drop=True)

df1 = df1[df1['Date'].diff().dt.total_seconds() > 0]

date_list = df1['Date'].values.astype('datetime64[s]').astype(float)

data = np.column_stack((df1['XGSM1'], df1['YGSM1'], df1['ZGSM1']))

store_data('X_Y_Z_gsm', data={'x': date_list, 'y': data})

pyspedas.kyoto.dst(trange=['1997-09-01', '1997-10-01'])
pyspedas.omni.data(trange=['1997-09-01', '1997-10-01'])

join_vec(['BX_GSE', 'BY_GSM', 'BZ_GSM'])

params = get_tsy_params('kyoto_dst',
'BX_GSE-BY_GSM-BZ_GSM_joined',
'proton_density',
'flow_speed',
't01', # or 't01', 'ts04'
pressure_tvar='Pressure',
speed=True)

tinterpol('X_Y_Z_gsm', 'proton_density')

tt01('X_Y_Z_gsm-itrp', parmod=params)

tplot('X_Y_Z_gsm-itrp_bt01')

The above code results in:

G variables required for T01 model; create a tplot variable containing the G variables, and provide the name of that keyword to the g_variables keyword.
tinterpol (linear) was applied to: X_Y_Z_gsm-itrp
parmod keyword required.
That name is currently not in pytplot
That name is currently not in pytplot
Variable not found: X_Y_Z_gsm-itrp_bt01

Thank you.

[jameswilburlewis]

I think the problem you're having with the T01 model is that you're not passing any values for the g1 and g2 model parameters it requires. (This is the g_variables argument to get_tsy_params)

The g1, g2 parameters for T01, and the w1...w6 parameters for TS04, can be calculated from OMNI data, but that code is not publicly available at the moment. (In IDL, we use Haje Korth's GEOPACK implementation that includes geopack_getg and geopack_getw routines which do exactly that, but that library is only distributed in binary form). The Python geopack package we're using in PySPEDAS does not offer these tools.

For TS04, it is possible to download yearly files with OMNI and w1...w6 coefficient values. PySPEDAS offers a get_w_params utility that can return the coefficient values using this source. There is also some Fortran code for generating these data files, which we can probably port to Python and add to the geopack library.

As far as I can tell, there's no such repository for the T01 g1, g2 coefficients, nor any Fortran code to calculate them. So for this model, we would either have to write a calculate_g_params routine, following the calculations outlined in the science publications describing the model, or possible get Haje to send us the Fortran code he's using for getpack_getg in the IDL GEOPACK library.

What we do in our test suite for T01 is to just supply some reasonable constants for g1 and g2 (just for the sake of exercising the code and verifying that IDL and Python are giving the same results). Here's what the test code looks like:

        gvars = np.zeros((len(mec.times), 2))
        # Note that mec.times may now have a different number of points than proton_density!
        gvars[:, 0] = np.repeat(6.0, len(mec.times))
        gvars[:, 1] = np.repeat(10.0, len(mec.times))
        store_data('g_variables', data={'x': mec.times, 'y': gvars})
        # This is a wrapper for get_tsy_params
        params = get_params('t01', g_variables='g_variables')

So we set g1 to 6.0 and g1 to 10.0, and just replicate those values to match the length of the input variables.

Maybe you can do something similar to at least get some test output from T01, even though the g1 and g2 values might not accurately reflect the actual conditions. In the meantime, we will try to implement equivalents to geopack_getg and geopack_getw, but that may take a bit of time.

I'm still looking at your CSV data file to see if I can figure out what's going on with the T96 model outputs.

[kmrs-iitd]

Thank you very much for taking your time and pointing out the issue. Actually, I'm currently studying magnetic field variations in the Earth's polar region using the T89, T96, T01, and TS04 models, and comparing them with observed data from the Polar satellite and SWMF output. If I use fixed g1=6 and g2=10, it could raise questions. Is the TS04 code also using arbitrary values for w1-w6 coefficients?

Another issue I encountered with the T96 and TS04 models is that I used these codes on one year of data averaged over 1 minute (i.e., 525,600 total data points), but I got 382,721 data points from TS04 and 377,492 from T96. Although, the input data contains no NaN or zero values.

[jameswilburlewis]

get_tsy_params will automatically download the w1...w6 parameters from Tsyganenko's site, so they should reflect the true conditions for the time range being loaded.

OMNI data can contain gaps, which would also cause gaps in the w1-w6 data. Perhaps that might explain some of the missing points. I can't say much more without knowing the specific times/positions/parameters for the missing output points....

[kmrs-iitd]

Thank you very much!
Someone recommended another method for calculating the external magnetic field using models like T89, T96, T01, TS05, and TS07. It seems to work with my dataset, but I'm not very familiar with it or its accuracy.
Link: https://spacepy.github.io/autosummary/spacepy.irbempy.get_Bfield.html#spacepy.irbempy.get_Bfield

I need assistance with another issue. I'm trying to increase the number of ticks on the y-axis using tplot_options('axis_tick_num', [(10,5)]), but it’s not having any effect.

Is there a way to plot multiple parameters on the x-axis similar to the example in the picture?

tplot_options('show_all_axes', True) is not working. Are there alternative ways to adjust y-axis ticks and plot multiple parameters on the x-axis, similar to the example in the attached picture?

pytplot.tplot("variable1", var_label = ["variable2", "variable3"]), this command gives error.

In given example it print few lines under x-axis with variable values.

Thank you

[jameswilburlewis]

Thanks very much for the link to the spacepy IRBEM routines! It might be very useful to provide a way for pyspedas to use that interface. It sounds like they are using compiled Fortran code, so it might be more efficient than the geopack package, at least for the models that IRBEM supports.

Would you mind filing your plotting questions as separate issues? It makes it easier to for us to track and assign tasks if it's just one topic per issue.

[kmrs-iitd]

Thanks , I will raise separate issues related to plotting.

[jameswilburlewis]

Were you able to resolve the issues you were having with unexpected model outputs and missing data points in the PySPEDAS geopack routines? Or did you decide to use the spacepy tools for your workflow?

[kmrs-iitd]

I am still unable to resolve the issues related to unexpected model outputs and missing data points in the PySPEDAS Geopack routines. I am using codes given in Geopack by Sheng Tian. In these codes (T89,T96, T01, T04), I am using Kp, Dst, Pres, IMF BY/Bz, G1, G2, and W1-W6 parameters from spacepy.omni.get_omni(df['Date'].tolist()), although it is taking much time to run.
One limitation of the SpacePy tools is that the code can only process 100,000 data points (d['ntime_max'] = 100000) . Please update when the issues with the PySPEDAS geopack routines have been resolved.