ancklo / ChaosMagPy

ChaosMagPy is a simple Python package for evaluating the CHAOS geomagnetic field model and other models of Earth's magnetic field.
Other
23 stars 3 forks source link

Broadcasting problem using synth_values #6

Closed sonjazw closed 3 years ago

sonjazw commented 3 years ago

Hi,

I would like to compute the field values for every point on a grid for a batch of coefficient. I can achieve it with a for loop and pass a single vector of coefficient and a meshgrid of theta and phi, but I was wondering if one could do it more elegantly by passing directly a multi-dimensional array of coefficients to "synth_values".

import numpy as np
from chaosmagpy.model_utils import synth_values

radius = 6371.2
theta = np.linspace(1., 179., num=179)
phi = np.linspace(-180., 180., num=361)

radius = radius*np.ones((phi.shape, ))

coeffs = np.random.rand(500, 195)

synth_values(coeffs, radius, theta, phi, nmax=13, mmax=13, grid=True) 

I keep getting the following error message :

ValueError: shape mismatch: objects cannot be broadcast to a single shape
ancklo commented 3 years ago

Hello,

if I understood correctly, you want to do something like the following?

radius = 6371.2
theta = np.linspace(1., 179., num=179)
phi = np.linspace(-180., 180., num=361)

# you don't need this, numpy will take care of it since a float can be broadcasted to any shape
# radius = radius*np.ones(phi.shape)

coeffs = np.random.rand(500, 1, 1, 195)  # <-- here I insert additional axis for numpy broadcasting

synth_values(coeffs, radius, theta, phi, nmax=13, mmax=13, grid=True)

This will return the vector components with shape (500, 179, 361). In other words, 500 maps of the magnetic field vector given on the theta-phi meshgrid.

You can also be more explicit by specifying the shapes of the arrays before passing them to synth_values (however now you should not use grid=True):

radius = 6371.2
theta = np.linspace(1., 179., num=179).reshape((1, 179, 1))
phi = np.linspace(-180., 180., num=361).reshape((1, 1, 361))
coeffs = np.random.rand(500, 1, 1, 195)

synth_values(coeffs, radius, theta, phi, nmax=13, mmax=13, grid=False)

However, I think the code becomes hard to read this way.

I hope this helps. Regards.

sonjazw commented 3 years ago

Hi,

Thank you very much for your answer ! It helped a lot !