"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Import notebook dependencies\n",
"\n",
"import sys\n",
"sys.path.append('..')\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn-white')\n",
"import numpy as np\n",
"import pandas as pd \n",
"from src import sha_lib as sha\n",
"import os"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. A more tangible example of modelling using spherical harmonics\n",
"\n",
"First, there's a data file with some spatial data to plot."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load the data\n",
"lsfile = '../data/external/land_5deg.csv'\n",
"lsdata = np.genfromtxt(lsfile, delimiter=',')\n",
"\n",
"x = lsdata[:,0].reshape(36, 72)\n",
"y = lsdata[:,1].reshape(36,72)\n",
"z = lsdata[:,2].reshape(36,72)\n",
"\n",
"# Plot them\n",
"plt.rcParams['figure.figsize'] = [15, 8]\n",
"plt.contourf(x,y,z)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data file consists of binary values (1=land, 0=water) on a 5-degree grid. We now take these data as the input to a global spherical harmonic analysis and calculate a spherical harmonic model with a user specified resolution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### >> USER INPUT HERE: Set the maximum spherical harmonic degree of the analysis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nmax = 5 # Max degree"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# These functions are used for the spherical harmonic analysis and synthesis steps\n",
"# Firstly compute a spherical harmonic model of degree and order nmax using the input data as plotted above\n",
"def sh_analysis(lsdata, nmax):\n",
" npar = (nmax+1)*(nmax+1)\n",
" ndat = len(lsdata)\n",
" lhs = np.zeros(npar*ndat).reshape(ndat,npar)\n",
"\n",
" rhs = np.zeros(ndat)\n",
" line = np.zeros(npar)\n",
" ic = -1\n",
" for i in range(ndat):\n",
" th = 90 - lsdata[i][1]\n",
" ph = lsdata[i][0] \n",
" rhs[i] = lsdata[i][2]\n",
" cmphi, smphi = sha.csmphi(nmax,ph)\n",
" pnm = sha.pnm_calc(nmax, th)\n",
" for n in range(nmax+1):\n",
" igx = sha.gnmindex(n,0)\n",
" ipx = sha.pnmindex(n,0)\n",
" line[igx] = pnm[ipx]\n",
" for m in range(1,n+1):\n",
" igx = sha.gnmindex(n,m)\n",
" ihx = sha.hnmindex(n,m)\n",
" ipx = sha.pnmindex(n,m)\n",
" line[igx] = pnm[ipx]*cmphi[m]\n",
" line[ihx] = pnm[ipx]*smphi[m]\n",
" lhs[i,:] = line\n",
"\n",
" shmod = np.linalg.lstsq(lhs, rhs.reshape(len(lsdata),1), rcond=None)\n",
" return(shmod)\n",
"\n",
"\n",
"# Now use the model to synthesise values on a 5 degree grid in latitude and longitude\n",
"def sh_synthesis(shcofs, nmax):\n",
" newdata =np.zeros(72*36*3).reshape(2592,3)\n",
" ic = 0\n",
" for ilat in range(36):\n",
" delta = 5*ilat+2.5\n",
" lat = 90 - delta\n",
" for iclt in range(72):\n",
" corr = 5*iclt+2.5\n",
" long = -180+corr\n",
" colat = 90-lat\n",
" cmphi, smphi = sha.csmphi(nmax,long)\n",
" vals = np.dot(sha.gh_phi(shcofs, nmax, cmphi, smphi), sha.pnm_calc(nmax, colat))\n",
" newdata[ic,0]=long\n",
" newdata[ic,1]=lat\n",
" newdata[ic,2]=vals\n",
" ic += 1\n",
" return(newdata)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Obtain the spherical harmonic coefficients\n",
"shmod = sh_analysis(lsdata=lsdata, nmax=nmax)\n",
"\n",
"# Read the model coefficients\n",
"shcofs = shmod[0]\n",
"\n",
"# Synthesise the model coefficients on a 5 degree grid\n",
"newdata = sh_synthesis(shcofs=shcofs, nmax=nmax)\n",
"# Reshape for plotting purposes\n",
"x = newdata[:,0].reshape(36, 72)\n",
"y = newdata[:,1].reshape(36,72)\n",
"z = newdata[:,2].reshape(36,72)\n",
"\n",
"# Plot the results\n",
"plt.rcParams['figure.figsize'] = [15, 8]\n",
"levels = [-1.5, -0.75, 0, 0.25, 0.5, 0.75, 2.]\n",
"plt.contourf(x,y,z)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2. Spherical harmonic models with data gaps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What happens if the data set is incomplete? In this section, you can experiment by removing data within a great circle of a specified radius and position. The functions below are used to create the data gap, and the modelling uses functions in the above section again."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def greatcircle(th1, ph1, th2, ph2):\n",
" th1 = np.deg2rad(th1)\n",
" th2 = np.deg2rad(th2)\n",
" dph = np.deg2rad(dlong(ph1,ph2))\n",
"\n",
" # Apply the cosine rule of spherical trigonometry\n",
" dth = np.arccos(np.cos(th1)*np.cos(th2) + \\\n",
" np.sin(th1)*np.sin(th2)*np.cos(dph))\n",
" return(dth)\n",
"\n",
"def dlong (ph1, ph2):\n",
" ph1 = np.sign(ph1)*abs(ph1)%360 # These lines return a number in the \n",
" ph2 = np.sign(ph2)*abs(ph2)%360 # range -360 to 360\n",
" if(ph1 < 0): ph1 = ph1 + 360 # Put the results in the range 0-360\n",
" if(ph2 < 0): ph2 = ph2 + 360\n",
" dph = max(ph1,ph2) - min(ph1,ph2) # So the answer is positive and in the\n",
" # range 0-360\n",
" if(dph > 180): dph = 360-dph # So the 'short route' is returned\n",
" return(dph)\n",
"\n",
"def gh_phi(gh, nmax, cp, sp):\n",
" rx = np.zeros(nmax*(nmax+3)//2+1)\n",
" igx=-1\n",
" igh=-1\n",
" for i in range(nmax+1):\n",
" igx += 1\n",
" igh += 1\n",
" rx[igx]= gh[igh]\n",
" for j in range(1,i+1):\n",
" igh += 2\n",
" igx += 1\n",
" rx[igx] = (gh[igh-1]*cp[j] + gh[igh]*sp[j])\n",
" return(rx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### >> USER INPUT HERE: Set the location and size of the data gap\n",
"\n",
"Colatitude: degrees\n",
"\n",
"Longitude: degrees\n",
"\n",
"Radius: km"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Remove a section of data centered on colat0, long0 and radius here\n",
"colat0 = 100\n",
"long0 = -55\n",
"radius = 500"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lsdata_gap = lsdata.copy()\n",
"for row in range(len(lsdata_gap)):\n",
" colat = 90 - lsdata_gap[row,1]\n",
" long = lsdata_gap[row,0]\n",
" if greatcircle(colat, long, colat0, long0) < radius/6371.2:\n",
" lsdata_gap[row,2] = np.nan\n",
"\n",
"print('Blanked out: ', np.count_nonzero(np.isnan(lsdata_gap)))\n",
"\n",
"x_gap = lsdata_gap[:,0].reshape(36, 72)\n",
"y_gap = lsdata_gap[:,1].reshape(36,72)\n",
"z_gap = lsdata_gap[:,2].reshape(36,72)\n",
"\n",
"# Plot the map with omitted data\n",
"plt.contourf(x_gap, y_gap, z_gap)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### >> USER INPUT HERE: Set the maximum spherical harmonic degree of the analysis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nmax = 5 #Max degree"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Select the non-nan data\n",
"lsdata_gap = lsdata_gap[~np.isnan(lsdata_gap[:,2])]\n",
"\n",
"# Obtain the spherical harmonic coefficients for the incomplete data set\n",
"shmod = sh_analysis(lsdata=lsdata_gap, nmax=nmax)\n",
"\n",
"# Read the model coefficients\n",
"shcofs = shmod[0]\n",
"\n",
"# Synthesise the model coefficients on a 5 degree grid\n",
"newdata = sh_synthesis(shcofs=shcofs, nmax=nmax)\n",
"# Reshape for plotting purposes\n",
"x_new = newdata[:,0].reshape(36, 72)\n",
"y_new = newdata[:,1].reshape(36,72)\n",
"z_new = newdata[:,2].reshape(36,72)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Print the maximum and minimum of the synthesised data. How do they compare to the original data, which was composed of only ones and zeroes? How do the max/min change as you vary the data gap size and the analysis resolution? Try this with a large gap, e.g. 5000 km."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"print(np.min(z_new))\n",
"print(np.max(z_new))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot the results with colour scale according to the data\n",
"plt.rcParams['figure.figsize'] = [15, 8]\n",
"plt.contourf(x_new, y_new, z_new)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now see what happens when we restrict the colour scale to values between 0 and 1 regardless of the data values (so that both ends of the colour scale are saturated)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot the results\n",
"plt.rcParams['figure.figsize'] = [15, 8]\n",
"levels = [0, 0.25, 0.5, 0.75, 1.]\n",
"plt.contourf(x_new, y_new, z_new, levels, extend='both')\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3. Making a spherical harmonic model of the geomagnetic field\n",
"\n",
"There is a file containing virtual observatory (VO) data calculated from the Swarm 3-satellite constellation in this repository. VOs use a local method to provide point estimates of the magnetic field at a specified time and fixed location at satellite altitude. This technique has various benefits, including: \n",
"\n",
"1. Easier comparison between satellite (moving instrument) and ground observatory (fixed instrument) data, which is particularly useful when studying time changes of the magnetic field, e.g. geomagnetic jerks.\n",
"2. One can compute VOs on a regular spatial grid to mitigate the effects of uneven ground observatory coverage.\n",
"\n",
"A brief summary of the method used to calculate each of these VO data points:\n",
"\n",
"1. Swarm track data over four months are divided into 300 globally distributed equal area bins\n",
"2. Gross outliers (deviating over 500nT from the CHAOS-6-x7 internal field model) are removed\n",
"3. Only data from magnetically quiet, dark times are kept\n",
"4. Estimates of the core and crustal fields from CHAOS-6-x7 are subtracted from each datum to give the residuals\n",
"5. A local magnetic potential $V$ is fit to the residuals within the cylinder\n",
"6. Values of the residual magnetic field in spherical coordinates ($r$, $\\theta$, $\\phi$) are computed from the obtained magnetic potential using $B=-\\nabla V$ at centre of the bin (at 490km altitude)\n",
"7.\tAn estimate of the modelled core field at the VO calculation point is added back onto the residual field to give the total internal magnetic field value\n",
"\n",
"In this section, we will use VOs as input data for a degree 13 spherical harmonic model of the geomagnetic field at 2015.0 and compare our (satellite data only) model to the latest IGRF (computed using both ground and satellite data) at the same time."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Read in the VO data\n",
"cfile = '../data/external/SwarmVO_IAGASummerSchool.dat'\n",
"cols = ['Year','Colat','Long','Rad','Br','Bt', 'Bp']\n",
"bvals = pd.read_csv(cfile, sep='\\s+', skiprows=0, header=None, index_col=0,\n",
" na_values=[99999.00,99999.00000], names=cols, comment='%')\n",
"bvals = bvals.dropna()\n",
"bvals['Bt']= -bvals['Bt']\n",
"bvals['Br']= -bvals['Br']\n",
"colnames=['Colat','Long','Rad','Bt','Bp', 'Br']\n",
"bvals=bvals.reindex(columns=colnames)\n",
"bvals.columns = ['Colat','Long','Rad','X','Y','Z']\n",
"\n",
"# Set the date to 2015.0 (this is the only common date between the VO data file and the IGRF coefficients file)\n",
"epoch = 2015.0\n",
"# Set the model resolution to 13 (the same as IGRF)\n",
"nmax = 13"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def geomagnetic_model(epoch, nmax, data, RREF=6371.2):\n",
" # Select the data for that date\n",
" colat = np.array(data.loc[epoch]['Colat'])\n",
" nx = len(colat) # Number of data triples\n",
" ndat = 3*nx # Total number of data\n",
" elong = np.array(data.loc[epoch]['Long'])\n",
" rad = np.array(data.loc[epoch]['Rad'])\n",
" rhs = data.loc[epoch][['X', 'Y', 'Z']].values.reshape(ndat,1)\n",
" npar = (nmax+1)*(nmax+1) # Number of model parameters\n",
" nx = len(colat) # Number of data triples\n",
" ndat = 3*nx # Total number of data\n",
"\n",
" # Arrays for the equations of condition\n",
" lhs = np.zeros(npar*ndat).reshape(ndat,npar)\n",
" linex = np.zeros(npar); liney = np.zeros(npar); linez = np.zeros(npar)\n",
"\n",
" iln = 0 # Row counter\n",
" for i in range(nx): # Loop over data triplets (X, Y, Z)\n",
" r = rad[i]\n",
" th = colat[i]; ph = elong[i] \n",
" rpow = sha.rad_powers(nmax, RREF, r)\n",
" cmphi, smphi = sha.csmphi(nmax, ph)\n",
" pnm, xnm, ynm, znm = sha.pxyznm_calc(nmax, th)\n",
"\n",
" for n in range(nmax+1): # m=0\n",
" igx = sha.gnmindex(n,0) # Index for g(n,0) in gh conventional order\n",
" ipx = sha.pnmindex(n,0) # Index for pnm(n,0) in array pnm\n",
" rfc = rpow[n]\n",
" linex[igx] = xnm[ipx]*rfc\n",
" liney[igx] = 0\n",
" linez[igx] = znm[ipx]*rfc\n",
"\n",
" for m in range(1,n+1): # m>0\n",
" igx = sha.gnmindex(n,m) # Index for g(n,m)\n",
" ihx = igx + 1 # Index for h(n,m)\n",
" ipx = sha.pnmindex(n,m) # Index for pnm(n,m)\n",
" cpr = cmphi[m]*rfc\n",
" spr = smphi[m]*rfc\n",
" linex[igx] = xnm[ipx]*cpr; linex[ihx] = xnm[ipx]*spr\n",
" liney[igx] = ynm[ipx]*spr; liney[ihx] = -ynm[ipx]*cpr\n",
" linez[igx] = znm[ipx]*cpr; linez[ihx] = znm[ipx]*spr\n",
"\n",
" lhs[iln, :] = linex\n",
" lhs[iln+1,:] = liney\n",
" lhs[iln+2 :] = linez\n",
" iln += 3\n",
"\n",
" shmod = np.linalg.lstsq(lhs, rhs, rcond=None) # Include the monopole\n",
" # shmod = np.linalg.lstsq(lhs[:,1:], rhs, rcond=None) # Exclude the monopole\n",
" gh = shmod[0]\n",
" return(gh)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"coeffs = geomagnetic_model(epoch=epoch, nmax=nmax, data=bvals)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Read in the IGRF coefficients for comparison"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"IGRF12_FILE = os.path.abspath('../data/external/igrf12coeffs.txt')\n",
"igrf12 = pd.read_csv(IGRF12_FILE, delim_whitespace=True, header=3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"igrf12.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Select the 2015 IGRF coefficients\n",
"igrf = np.array(igrf12[str(epoch)])\n",
"# Insert a zero in place of the monopole coefficient (this term is NOT included in the IRGF, or any other field model)\n",
"igrf = np.insert(igrf, 0, 0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compare our model coefficients to IGRF-12"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"print(\"VO only model \\t IGRF12\")\n",
"for i in range(196):\n",
" print(\"% 9.1f\\t%9.1f\" %(coeffs[i], igrf[i]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"The only common date between the VO data file and the IGRF coefficients file is 2015.0. Use VOs to produce a model at a time of your choice between 2014 and 2018 (other than 2015.0 as we did above) and compare your coefficients to the IGRF at that time. **Hints:** For dates other than yyyy.0, you will need to interpolate the VO data. You may also need to interpolate the IGRF coefficients between 2010.0 and 2015.0. The IGRF coefficients between 2015.0 and 2020.0 will need to be extrapolated from the 2015.0 values using the given secular variation predictions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"Create a second spherical harmonic model based on ground observatory data only and compare it to your satellite data only model and to the IGRF. To do this, you can use observatory annual means as the input data, which we have supplied in the file `oamjumpsapplied.dat` in the external data directory. **Hint:** The `mag` module contains a function that will extract all annual mean values for all observatories in a given year. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Acknowledgements\n",
"\n",
"The land/sea data file was provided by John Stevenson (BGS).\n",
"The VO data were computed by Magnus Hammer (DTU) and the file prepared by Will Brown (BGS)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}