Internal vs External SEC Geometry
Closed this issue · 10 comments
I've noticed some unexpected behavior while trying to decompose ground magnetic field observations into a combination of external and internal SECS (as in the Pulkkinen et al. EPS 2003 paper). With two sets of SECS, one external shell (at say R_Earth + 100km) and one internal shell (at say R_Earth - 30km), I would expect the resulting current vectors to be generally pointing in the opposite direction (since currents have to move in the opposite direction across the observation plane in order to produce magnetic fields pointing in the same direction on that plane). However, the SECS decomposition algorithm returns currents that point in the same direction regardless of the shell height with respect to observations (so regardless of whether the shell is positioned above or below the Earth's surface).
This synthetic example (adapted from the notebook example) demonstrates the problem:
import numpy as np
from pySECS import SECS
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
R_E = 6371.e3
sec_height = -30.e3
# place current system under Golden CO
sec_loc = np.array([39.7555, -105.2211, R_E + sec_height])
# set up divergence-free current system
system_df = SECS(sec_df_loc=sec_loc)
system_df.fit_unit_currents()
system_df.sec_amps *= 10000.
# prediction grid
lats = np.arange(30.5, 50.5, 1.)
lons = np.arange(-125.5, -85.5, 1.)
lat_secs, lon_secs = np.meshgrid(lats, lons)
lat_secs = lat_secs.flatten()
lon_secs = lon_secs.flatten()
secs_current_coords = np.column_stack((lat_secs, lon_secs, (R_E + sec_height) * np.ones(lat_secs.size)))
secs_bfield_coords = np.column_stack((lat_secs, lon_secs, R_E * np.ones(lat_secs.size)))
# predict stuff
j = system_df.predict_J(secs_current_coords)
b = system_df.predict_B(secs_bfield_coords)
# plot
scale = '10m'
coast = cfeature.NaturalEarthFeature(category='physical', scale=scale, edgecolor='k',
facecolor='none', name='coastline')
countries = cfeature.NaturalEarthFeature(category='cultural', scale=scale, edgecolor='k',
facecolor='none', name='admin_0_countries')
projection = ccrs.LambertConformal(central_longitude=-96.0, central_latitude=39.0,
standard_parallels=(33., 45.))
proj_data = ccrs.PlateCarree()
fig = plt.figure(figsize=(20., 14.))
ax_j = fig.add_subplot(121, projection=projection)
ax_j.set_title("Currents @ {:.0f} km".format(sec_height/1.e3))
ax_j.add_feature(coast)
ax_j.add_feature(countries)
ax_j.set_xlim([-3.e6, 3.e6])
ax_j.set_ylim([-2.e6, 2.e6])
ax_j.quiver(secs_current_coords[:, 1], secs_current_coords[:, 0], j[:, 1], j[:, 0],
transform=proj_data, regrid_shape=30, color='r', zorder=9)
ax_b = fig.add_subplot(122, projection=projection)
ax_b.set_title("Surface Magnetic Fields")
ax_b.add_feature(coast)
ax_b.add_feature(countries)
ax_b.set_xlim([-3.e6, 3.e6])
ax_b.set_ylim([-2.e6, 2.e6])
ax_b.quiver(secs_bfield_coords[:, 1], secs_bfield_coords[:, 0], b[:, 1], b[:, 0],
transform=proj_data, regrid_shape=30, color='b', zorder=9)
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=0.95)
plt.show()The resulting output looks like this:
For counter-clockwise current loops underneath the observer (so within the Earth, r < R_Earth), the magnetic fields should be pointing radially outward, the opposite of what the algorithm returns. The output at the Earth's surface makes sense when SECS are placed outside the Earth (r > R_Earth), as the predicted magnetic fields point radially inwards. But B still points radially inwards for SECS placed inside the Earth. In general, for this single SEC example, B points radially inwards regardless of where the SEC is placed with respect to the observation shell.
I think the problem might have to do with the last few lines of the T_df function in secs.py (lines 295-297):
https://github.com/greglucas/pySECS/blob/198d8151791c9910c302b1fbba85eb61bf971c87/pySECS/secs.py#L292-L297
Since the equations for B_theta are in a spherical coordinate system with respect to the SEC pole, you have to do something to get B_theta into the observation (global, geomagnetic) coordinate system, hence these few lines of code. It looks like the factor of -1 in lines 295 and 296 (on the Bx and By components) are intended to flip the decomposed B_theta to point towards the SEC pole. This will be correct for SECS above observation locations (e.g., above the Earth's surface) but will not be correct for SECS below the observation locations.
I might be misunderstanding what's going on in the code and therefore be wrong about this, but I do get results that look correct (with the radial direction of B changing when moving the SECS from outside to inside the Earth) when I change this geometry factor in my copy of the code. I think the change at the end of the T_df function should be fairly simple:
geometry_factor = - np.ones(Btheta.shape)
under_locs = sec_r < obs_r
geometry_factor[under_locs] = 1.
# Transform back to Bx, By, Bz at each local point
T = np.empty((nobs, 3, nsec))
# alpha == angle (from cartesian x-axis (By), going towards y-axis (Bx))
T[:, 0, :] = geometry_factor * Btheta * np.sin(alpha)
T[:, 1, :] = geometry_factor * Btheta * np.cos(alpha)
T[:, 2, :] = -Br(Br should still be multiplied by -1 to swap between r being radially outward from the Earth and geomagnetic Z pointing towards the center of the Earth.)
If you think I'm correct about all this, happy to take out a PR with the fix. And thanks for nicely coding up this algorithm -- it's been wonderfully quick and easy to get this working.
@bmurphy-usgs, that is an excellent write-up! Admittedly, I have not used interior SECS, so I did not pay as much attention to that use case (other than putting the flag in). I think you are correct that due to my "simplifying" angle transforms with the directions that the factors should change when going through the plane of data.
It would be most excellent if you could make a PR with this change. Note that @erigler-usgs might be interested in this update as well and have comments about it. I know he was in the process of making an official release of pySECS through the USGS.
I would only suggest that, prior to accepting @bmurphy-usgs's PR, we test against my original code at https://github.com/usgs/geomag-imp. Not because I trust my own ability to do spherical geometry properly, but because I did at least test against Antti Pulkkinnen's original Matlab code, which I'm pretty sure was once checked thoroughly by Olaf Amm, the original author of this technique. I'm happy to help make this happen.
@erigler-usgs, that sounds good, I'll work on getting your original version working with my test cases and let you know if I have any questions
I just went through and checked the equations from the original paper
a little more thoroughly and if you look at equation A.8, I am missing the minus sign in that calculation:
https://github.com/greglucas/pySECS/blob/198d8151791c9910c302b1fbba85eb61bf971c87/pySECS/secs.py#L284
So, rather than needing a geometry factor, I think the negative just needs to be put on that line. Ben could you test that fix out and make that small PR if that solves the issue? Then we can address other tests of codes into a separate PR if more errors are found. But, after reading through the paper again I'm pretty sure I just mistyped that equation...
Sure enough, that fix will probably have the same effect as my hack-y -1 multiplication. Thought I had checked equations A.7 and A.8, but I guess I just didn't pay close enough attention. I'll try this and see how things look.
Yep that seems to have fixed it. And a better fix than my hacky work around. I'll take out a PR to update this. Should I also update the version tag to 0.0.2?
Yes, please feel free to update the version tag. Then once we get some more concrete tests added in we can make a 0.1 release.
@bmurphy-usgs I just pushed a lot of new stuff up so you will need to rebase when making the PR. I added a test that is currently failing and will hopefully pass with your change to close this issue.
@bmurphy-usgs just pinging you to see if you've had a chance to look into the PR any more?
@greglucas, quick PR (#3) fixes the missing minus sign issue, which does seem to fix this internal vs external problem. If the PR looks OK, we can close this issue. But I'd also like to get in some more rigorous tests comparing @erigler-usgs's original SECS code to this code sometime soon under an additional PR. Haven't really looked at the tests you've implemented yet, but this could go along with those (I'm envisioning some direct comparison of output between the two codes).