Notebook

Test of Magnetic Model Evaluation - CHAOS-6-Core¶

In [1]:

from time import time as get_time

class Timer():
    def __init__(self):
        self.start = get_time()
 
    def __call__(self, label="elapsed time", restart=False):
        print("%s %gs" % (label, get_time() - self.start))
        if restart:
            self.start = get_time()

Load model¶

In [2]:

import sys
from os.path import basename
from eoxmagmod.data import CHAOS6_CORE_LATEST as MODEL_FILE

print(basename(MODEL_FILE))

with open(MODEL_FILE) as fin:
    for idx, line in enumerate(fin, 1):
        if idx > 5:
            print("...")
            break
        sys.stdout.write(line)

CHAOS-6-x7_core.shc
# CHAOS-6-x7 core field model                      
# Based on model iteration_1 of CF_CHAOS_09d_18.mat
# extracted on 09-Oct-2018 12:11:25                
  1  20 221   6   5
          1997.1020   1997.2019   1997.3018   1997.4018   1997.5017   1997.6016   1997.7016   1997.8015   1997.9014   1998.0014   1998.1013   1998.2012   1998.3012   1998.4011   1998.5010   1998.6010   1998.7009   1998.8008   1998.9008   1999.0007   1999.1006   1999.2005   1999.3005   1999.4004   1999.5003   1999.6003   1999.7008   1999.8012   1999.9017   2000.0022   2000.1027   2000.2026   2000.3025   2000.4025   2000.5024   2000.6023   2000.7023   2000.8022   2000.9021   2001.0021   2001.1020   2001.2019   2001.3018   2001.4018   2001.5017   2001.6016   2001.7016   2001.8015   2001.9014   2002.0014   2002.1013   2002.2012   2002.3012   2002.4011   2002.5010   2002.6010   2002.7009   2002.8008   2002.9008   2003.0007   2003.1006   2003.2005   2003.3005   2003.4004   2003.5003   2003.6003   2003.7008   2003.8012   2003.9017   2004.0022   2004.1027   2004.2026   2004.3025   2004.4025   2004.5024   2004.6023   2004.7023   2004.8022   2004.9021   2005.0021   2005.1020   2005.2019   2005.3018   2005.4018   2005.5017   2005.6016   2005.7016   2005.8015   2005.9014   2006.0014   2006.1013   2006.2012   2006.3012   2006.4011   2006.5010   2006.6010   2006.7009   2006.8008   2006.9008   2007.0007   2007.1006   2007.2005   2007.3005   2007.4004   2007.5003   2007.6003   2007.7008   2007.8012   2007.9017   2008.0022   2008.1027   2008.2026   2008.3025   2008.4025   2008.5024   2008.6023   2008.7023   2008.8022   2008.9021   2009.0021   2009.1020   2009.2019   2009.3018   2009.4018   2009.5017   2009.6016   2009.7016   2009.8015   2009.9014   2010.0014   2010.1013   2010.2012   2010.3012   2010.4011   2010.5010   2010.6010   2010.7009   2010.8008   2010.9008   2011.0007   2011.1006   2011.2005   2011.3005   2011.4004   2011.5003   2011.6003   2011.7008   2011.8012   2011.9017   2012.0022   2012.1027   2012.2026   2012.3025   2012.4025   2012.5024   2012.6023   2012.7023   2012.8022   2012.9021   2013.0021   2013.1020   2013.2019   2013.3018   2013.4018   2013.5017   2013.6016   2013.7016   2013.8015   2013.9014   2014.0014   2014.1013   2014.2012   2014.3012   2014.4011   2014.5010   2014.6010   2014.7009   2014.8008   2014.9008   2015.0007   2015.1006   2015.2005   2015.3005   2015.4004   2015.5003   2015.6003   2015.7008   2015.8012   2015.9017   2016.0022   2016.1027   2016.2026   2016.3025   2016.4025   2016.5024   2016.6023   2016.7023   2016.8022   2016.9021   2017.0021   2017.1020   2017.2019   2017.3018   2017.4018   2017.5017   2017.6016   2017.7016   2017.8015   2017.9014   2018.0014   2018.1013   2018.2012   2018.3012   2018.4011   2018.5010   2018.6010   2018.7009   2018.8008   2018.9008   2019.0007   2019.1006
...

/home/pacesm/conda/lib/python3.7/site-packages/spacepy/pycdf/__init__.py:1209: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working
  class CDF(collections.MutableMapping):

In [3]:

# generate random points
from math import pi
from numpy import sin, cos, dstack, linspace
from numpy.random import random, uniform

EARTH_RADIUS = 6371.2 # km
N_coords = 500

X = random((N_coords*2, 2))
X = X[X[..., 1] < sin(X[..., 0]*pi)][:N_coords, ...]
X[..., 0] = pi * (X[..., 0] - 0.5)
X[..., 1] = pi*(2.0*X[..., 1] / cos(X[..., 0]) - 1.0)
X *= 180./pi
lats = X[..., 0]
lons = X[..., 1]
rads = uniform(EARTH_RADIUS, 1.1*EARTH_RADIUS, N_coords)

coords = dstack((lats, lons, rads))[0]

In [4]:

# times
from numpy import linspace
from eoxmagmod import load_model_shc

model = load_model_shc(MODEL_FILE)
start, end = model.validity

N_times = 1500

times = linspace(start, end, N_times)

default version¶

interpolated in the MJD2000 time domain
decimal years aligned with callendar years

In [5]:

from numpy import empty
from eoxmagmod import load_model_shc

model_mjd2k = load_model_shc(MODEL_FILE)
print(model_mjd2k.coefficients)


timer = Timer()

b_nec_mjd2k = empty((N_times, N_coords, 3))
for idx, time in enumerate(times):
    b_nec_mjd2k[idx, ...] = model_mjd2k.eval(time, coords, scale=[1, 1, -1])

timer("model evaluated in")

<eoxmagmod.magnetic_model.coefficients.SparseSHCoefficientsTimeDependent object at 0x7fb47472f780>
model evaluated in 1.67441s

decimal year interpolation¶

intrepolated in the decimal year time domain
decimal years aligned with callendar years

In [6]:

from numpy import empty
from eoxmagmod import load_model_shc

model_dy = load_model_shc(
    MODEL_FILE, interpolate_in_decimal_years=True
)
print(model_dy.coefficients)

timer = Timer()

b_nec_dy = empty((N_times, N_coords, 3))
for idx, time in enumerate(times):
    b_nec_dy[idx, ...] = model_dy.eval(time, coords, scale=[1, 1, -1])

timer("model evaluated in")

<eoxmagmod.magnetic_model.coefficients.SparseSHCoefficientsTimeDependentDecimalYear object at 0x7fb44d1430f0>
model evaluated in 1.68006s

simplified MJD2000 calculation¶

intrepolated in the MJD2000 time domain
decimal years calculated using 365.25 days per year
this method is used by the chaosmagpy

In [7]:

from numpy import empty
from eoxmagmod import load_model_shc
from eoxmagmod.time_util import decimal_year_to_mjd2000_simple

model_mjd2k_simple = load_model_shc(
    MODEL_FILE, to_mjd2000=decimal_year_to_mjd2000_simple,
)
print(model_mjd2k_simple.coefficients)

timer = Timer()

b_nec_mjd2k_simple = empty((N_times, N_coords, 3))
for idx, time in enumerate(times):
    b_nec_mjd2k_simple[idx, ...] = model_mjd2k_simple.eval(time, coords, scale=[1, 1, -1])

timer("model evaluated in")

<eoxmagmod.magnetic_model.coefficients.SparseSHCoefficientsTimeDependent object at 0x7fb44d143438>
model evaluated in 1.67169s

Plotting the results¶

In [8]:

from numpy import zeros, concatenate, isnan
from matplotlib import pyplot as plt
from eoxmagmod import vnorm
from eoxmagmod.time_util import mjd2000_to_decimal_year_simple as mjd2000_to_decimal_year


def plot(ax, x0, y0, label):
    x, y = [], []
    for idx in range(0, N_coords):
        x.append(x0)
        y.append(y0[:, idx])
    x, y = concatenate(x), concatenate(y)
    plt.plot(x, y, '.', markersize=2, color='#1f77b4', label='random samples')
    plt.plot(node_times, node_val, '+r', markersize=5, label='interpolation nodes')
    plt.ylabel(label)
    plt.legend(loc='upper right')
    plt.text(
        0.01, 0.95, 'max: %g nT' % y[~isnan(y)].max(),
        horizontalalignment='left', verticalalignment='top', transform=ax.transAxes
    )
    plt.text(
        0.01, 0.05, 'min: %g nT' % y[~isnan(y)].min(),
        horizontalalignment='left', verticalalignment='bottom', transform=ax.transAxes
    )
    plt.grid()

Effect of the MJD2k vs. decimal_year interpolation time domain¶

In [9]:

model = load_model_shc(MODEL_FILE)
node_times = mjd2000_to_decimal_year(model.coefficients._times)
node_val = zeros(node_times.shape)

delta_b_nec = b_nec_dy - b_nec_mjd2k
delta_f = vnorm(b_nec_dy) - vnorm(b_nec_mjd2k)

fig = plt.figure(figsize=(16, 16), dpi=100)
plt.title("decimal year vs. MJD2000 time domain")

times_dy = mjd2000_to_decimal_year(times)
    
ax = plt.subplot(4, 1, 1)
plot(ax, times_dy, delta_b_nec[..., 0], "delta B_N / nT")

ax = plt.subplot(4, 1, 2)
plot(ax, times_dy, delta_b_nec[..., 1], "delta B_E / nT")

ax = plt.subplot(4, 1, 3)
plot(ax, times_dy, delta_b_nec[..., 2], "delta B_C / nT")

ax = plt.subplot(4, 1, 4)
plot(ax, times_dy, delta_f, "delta F / nT")

plt.show()
%matplotlib inline

Effect of the regular vs. simplified decimal year calculation¶

In [10]:

model = load_model_shc(MODEL_FILE)
node_times = mjd2000_to_decimal_year(model.coefficients._times)
node_val = zeros(node_times.shape)

delta_b_nec = b_nec_mjd2k_simple - b_nec_mjd2k
delta_f = vnorm(b_nec_mjd2k_simple) - vnorm(b_nec_mjd2k)

fig = plt.figure(figsize=(16, 16), dpi=100)
plt.title("decimal year vs. MJD2000 time domain")

times_dy = mjd2000_to_decimal_year(times)

ax = plt.subplot(4, 1, 1)
plot(ax, times_dy, delta_b_nec[..., 0], "delta B_N / nT")

ax = plt.subplot(4, 1, 2)
plot(ax, times_dy, delta_b_nec[..., 1], "delta B_E / nT")

ax = plt.subplot(4, 1, 3)
plot(ax, times_dy, delta_b_nec[..., 2], "delta B_C / nT")

ax = plt.subplot(4, 1, 4)
plot(ax, times_dy, delta_f, "delta F / nT")

plt.show()
%matplotlib inline

Plot the random locations¶

In [11]:

from cartopy.feature import LAND, OCEAN, COASTLINE
from cartopy.crs import Mollweide, PlateCarree
from matplotlib import pyplot as plt
from matplotlib.cm import winter as colormap
from matplotlib.colors import Normalize
from matplotlib.colorbar import ColorbarBase
from mpl_toolkits.mplot3d import Axes3D
from eoxmagmod import convert, GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN
%matplotlib inline

norm = Normalize(vmin=1,vmax=1.1)

#help(ccrs)
fig = plt.figure(figsize=(16, 16), dpi=100, facecolor='w', edgecolor='k')
ax = plt.subplot(2, 1, 1, projection=Mollweide())
gl = ax.gridlines(crs=PlateCarree(), draw_labels=False, linewidth=1, color='silver', alpha=0.5, linestyle='--')
ax.add_feature(LAND, facecolor=(1.0, 1.0, 0.9))
ax.add_feature(OCEAN, facecolor=(0.9, 1.0, 1.0))
ax.add_feature(COASTLINE, edgecolor='silver')


obj = ax.scatter(
    lons, lats, c=rads/EARTH_RADIUS, s=1.5,
    cmap=colormap, norm=norm, transform=PlateCarree(),
     
)
cb = plt.colorbar(obj)
cb.ax.set_ylabel("r/R (R = %g km)" % EARTH_RADIUS)


ax = plt.subplot(2, 1, 2, projection='3d')

cart_coords = convert(coords, GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN)/EARTH_RADIUS
obj = ax.scatter(
    cart_coords[:, 0], cart_coords[:, 1], cart_coords[:, 2], 
    c=rads/EARTH_RADIUS, s=1.5,
    cmap=colormap, norm=norm,
)

ax.set_aspect('equal','box')
ax.set_xlabel("x/R")
ax.set_ylabel("y/R")
ax.set_zlabel("z/R")

plt.show()