Sampling data from a 3D model

In this example we’ll see how to sample a 3D model output at arbitrary points within the model domain.

First, load the required modules.

from psipy.model import MASOutput
from psipy.data import sample_data

import astropy.constants as const
import astropy.units as u

import matplotlib.pyplot as plt
import numpy as np

Load a set of MAS output files, and get the number density variable from the model run.

mas_path = sample_data.mas_helio()
model = MASOutput(mas_path)
rho = model['rho']

Out:

Files Downloaded:   0%|          | 0/3 [00:00<?, ?file/s]

rho002.hdf:   0%|          | 0.00/8.02M [00:00<?, ?B/s]


vr002.hdf:   0%|          | 0.00/7.96M [00:00<?, ?B/s]



br002.hdf:   0%|          | 0.00/7.94M [00:00<?, ?B/s]

rho002.hdf:   0%|          | 100/8.02M [00:00<2:48:31, 793B/s]


vr002.hdf:   0%|          | 100/7.96M [00:00<2:52:23, 769B/s]

rho002.hdf:   1%|          | 51.7k/8.02M [00:00<00:32, 248kB/s]


vr002.hdf:   0%|          | 11.2k/7.96M [00:00<02:34, 51.5kB/s]

rho002.hdf:   3%|3         | 266k/8.02M [00:00<00:08, 903kB/s]


vr002.hdf:   1%|          | 53.2k/7.96M [00:00<00:45, 175kB/s]

rho002.hdf:   9%|9         | 761k/8.02M [00:00<00:03, 2.12MB/s]


vr002.hdf:   2%|1         | 144k/7.96M [00:00<00:20, 384kB/s]

rho002.hdf:  21%|##1       | 1.69M/8.02M [00:00<00:01, 4.32MB/s]


vr002.hdf:   5%|4         | 385k/7.96M [00:00<00:08, 913kB/s]

rho002.hdf:  33%|###3      | 2.66M/8.02M [00:00<00:00, 5.94MB/s]

rho002.hdf:  43%|####3     | 3.45M/8.02M [00:00<00:00, 5.94MB/s]


vr002.hdf:   8%|8         | 656k/7.96M [00:00<00:06, 1.11MB/s]

rho002.hdf:  51%|#####     | 4.07M/8.02M [00:01<00:01, 3.69MB/s]


vr002.hdf:  10%|9         | 762k/7.96M [00:01<00:11, 648kB/s]

rho002.hdf:  59%|#####9    | 4.73M/8.02M [00:01<00:00, 4.25MB/s]


vr002.hdf:  15%|#4        | 1.16M/7.96M [00:01<00:05, 1.20MB/s]

rho002.hdf:  66%|######5   | 5.26M/8.02M [00:01<00:00, 4.21MB/s]


vr002.hdf:  17%|#6        | 1.34M/7.96M [00:01<00:05, 1.23MB/s]

rho002.hdf:  72%|#######1  | 5.76M/8.02M [00:01<00:00, 4.08MB/s]


vr002.hdf:  20%|#9        | 1.57M/7.96M [00:01<00:04, 1.46MB/s]

rho002.hdf:  78%|#######7  | 6.22M/8.02M [00:01<00:00, 3.97MB/s]

rho002.hdf:  83%|########2 | 6.65M/8.02M [00:01<00:00, 3.36MB/s]

rho002.hdf:  88%|########7 | 7.02M/8.02M [00:01<00:00, 3.12MB/s]

rho002.hdf:  92%|#########1| 7.36M/8.02M [00:02<00:00, 2.43MB/s]



br002.hdf:   0%|          | 100/7.94M [00:02<50:33:19, 43.6B/s]

rho002.hdf:  95%|#########5| 7.64M/8.02M [00:02<00:00, 2.41MB/s]



br002.hdf:   0%|          | 10.7k/7.94M [00:02<21:18, 6.20kB/s]



br002.hdf:   1%|          | 51.3k/7.94M [00:02<03:34, 36.7kB/s]



br002.hdf:   2%|2         | 164k/7.94M [00:02<00:55, 141kB/s]


vr002.hdf:  22%|##2       | 1.76M/7.96M [00:02<00:13, 476kB/s]



br002.hdf:   7%|6         | 535k/7.94M [00:02<00:12, 572kB/s]


vr002.hdf:  24%|##4       | 1.92M/7.96M [00:02<00:10, 556kB/s]



br002.hdf:  14%|#4        | 1.13M/7.94M [00:02<00:05, 1.35MB/s]


vr002.hdf:  29%|##8       | 2.29M/7.96M [00:02<00:06, 899kB/s]


vr002.hdf:  37%|###6      | 2.92M/7.96M [00:03<00:03, 1.60MB/s]

rho002.hdf:  99%|#########8| 7.90M/8.02M [00:03<00:00, 1.05MB/s]


Files Downloaded:  33%|###3      | 1/3 [00:03<00:07,  3.51s/file]



br002.hdf:  20%|##        | 1.59M/7.94M [00:03<00:06, 1.02MB/s]


vr002.hdf:  40%|####      | 3.22M/7.96M [00:03<00:04, 1.10MB/s]


vr002.hdf:  43%|####3     | 3.45M/7.96M [00:03<00:03, 1.16MB/s]



br002.hdf:  23%|##2       | 1.82M/7.94M [00:03<00:06, 944kB/s]


vr002.hdf:  48%|####7     | 3.80M/7.96M [00:03<00:02, 1.44MB/s]



br002.hdf:  27%|##7       | 2.16M/7.94M [00:03<00:04, 1.21MB/s]


vr002.hdf:  57%|#####6    | 4.51M/7.96M [00:03<00:01, 2.35MB/s]



br002.hdf:  35%|###4      | 2.74M/7.94M [00:03<00:02, 1.87MB/s]


vr002.hdf:  64%|######3   | 5.07M/7.96M [00:04<00:00, 2.95MB/s]



br002.hdf:  42%|####2     | 3.34M/7.94M [00:04<00:01, 2.57MB/s]


vr002.hdf:  74%|#######4  | 5.91M/7.96M [00:04<00:00, 4.01MB/s]



br002.hdf:  47%|####7     | 3.75M/7.94M [00:04<00:01, 2.87MB/s]


vr002.hdf:  87%|########7 | 6.93M/7.96M [00:04<00:00, 5.41MB/s]



br002.hdf:  56%|#####6    | 4.48M/7.94M [00:04<00:00, 3.73MB/s]



Files Downloaded:  67%|######6   | 2/3 [00:04<00:02,  2.01s/file]



br002.hdf:  62%|######2   | 4.96M/7.94M [00:04<00:01, 1.84MB/s]



br002.hdf:  67%|######7   | 5.33M/7.94M [00:05<00:01, 1.97MB/s]



br002.hdf:  74%|#######4  | 5.88M/7.94M [00:05<00:00, 2.50MB/s]



br002.hdf:  85%|########4 | 6.75M/7.94M [00:05<00:00, 3.60MB/s]



br002.hdf:  92%|#########2| 7.33M/7.94M [00:05<00:00, 3.89MB/s]



br002.hdf:  99%|#########8| 7.85M/7.94M [00:05<00:00, 3.02MB/s]




Files Downloaded: 100%|##########| 3/3 [00:05<00:00,  1.69s/file]
Files Downloaded: 100%|##########| 3/3 [00:05<00:00,  1.92s/file]

Choose a set of 1D points to interpolate the model output at.

Here we keep a constant radius, and a set of longitudes that go all the way from 0 to 360 degrees. Then we choose two different, but close latitude values, and plot the results.

As expected, the values at 0 and 360 degrees are the same, and the profiles are similar, but different, due to the small difference in latitude.

fig, ax = plt.subplots()

npoints = 1000
r = 50 * np.ones(npoints) * const.R_sun
lon = np.linspace(0, 360, npoints) * u.deg

for latitude in [0, 1] * u.deg:
    lat = latitude * np.ones(npoints)
    samples = rho.sample_at_coords(lon, lat, r)

    ax.plot(lon, samples, label='lat = ' + str(latitude))

ax.legend()
ax.set_xlim(0, 360)
ax.set_ylim(bottom=0)
ax.set_xlabel('Longitude (deg)')
ax.set_ylabel(r'$\rho$ (cm$^{-3}$)')
ax.set_xticks([0, 90, 180, 270, 360])
plt.show()