{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# IBP-Model\n", "\n", "> Notebook authors: Ina Rusch\n", "\n", "> Package authors: Martin Rother, Ina Rusch" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Resources\n", "\n", "Python package: `ibpmodel` \n", "[Code](https://igit.iap-kborn.de/ibp/ibp-model) \n", "[Documentation](https://ibp-model.readthedocs.io/) \n", "[Swarm Data Handbook](https://swarmhandbook.earth.esa.int/catalogue/SW_IBP_CLI_2_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic usage" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import ibpmodel as ibp\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the IBP model at given time of year, position, and conditions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.calculateIBPindex(\n", " day_month=15, # Day of Year or Month\n", " longitude=0, # Longitude in degree\n", " local_time=20.9, # Local time in hours\n", " f107=150, # F10.7 cm Solar Flux index\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.calculateIBPindex(day_month=['Jan','Feb','Mar'], local_time=22)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the global pattern for a given day:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.plotIBPindex(doy=349)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the month vs longitude pattern:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.plotButterflyData(f107=150)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More examples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.calculateIBPindex('Jan', longitude=170, local_time=np.arange(-2,-1,0.1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.calculateIBPindex(day_month=[1,15,31], longitude=[-170,175,170], local_time=0, f107=150)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.plotIBPindex('Feb', cmap='viridis', alpha=0.9)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "doys = [319, 349, 15, 46]\n", "fig_doys, axes_doys = plt.subplots(len(doys),1, layout='constrained',figsize=(9, 15))\n", " \n", "for d, ax in zip(doys, axes_doys):\n", " ax, scalarmap = ibp.plotIBPindex(d, ax=ax)\n", "\n", "ibp.ibpforward.setcolorbar(scalarmap, fig_doys, axes_doys, fraction=0.05)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solar = [90, 130, 180]\n", "doy = 349\n", "fig_solar, axes_solar = plt.subplots(len(solar),1, layout='constrained',figsize=(9, 13))\n", "fig_solar.suptitle(f\"Doy {doy} with different solar radio flux\", fontsize=16)\n", "\n", "for f, ax in zip(solar, axes_solar):\n", " ax, scalarmap = ibp.plotIBPindex(doy, f107=f, ax=ax)\n", " ax.set_title(f\"IBP by solar radio flux {f}\")\n", "\n", "cbar = ibp.ibpforward.setcolorbar(scalarmap, fig_solar, axes_solar, fraction=0.05, orientation='horizontal')\n", "cbar.set_label('new label name')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ibp.plotButterflyData(150, cmap='magma')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solar = [130, 180]\n", "fig_bfly, axes_bfly = plt.subplots(1,len(solar), layout='constrained',figsize=(12.4, 5))\n", " \n", "for b, ax in zip(solar, axes_bfly):\n", " ax, scalarmap = ibp.plotButterflyData(b, ax=ax)\n", "\n", "ibp.ibpforward.setcolorbar(scalarmap, fig_bfly, axes_bfly)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualisation along Swarm orbits\n", "\n", "Source: \n", "https://igit.iap-kborn.de/ibp/ibp-model/-/blob/main/template/IBP-VirES.ipynb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install --quiet lxml" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import ibpmodel as ibp\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from datetime import datetime, timezone, timedelta, date\n", "from viresclient import SwarmRequest" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def getlt(datum, lon):\n", " '''Output of the day of the year and local time for specifying date(s) and longitude(s).\n", " \n", " Parameters\n", " ----------\n", " datum : datetime or list of datetimes\n", " lon : int or list of int\n", " The geographical longitude, ``-180 <= longitude <= 180``.\n", " \n", " Returns\n", " -------\n", " tupel\n", " contains (list of) day(s) of the year and (list of) local time(s)\n", " '''\n", " def calclt(date, lon):\n", " return round((date.hour + date.minute/60 + date.second/3600 + 24*lon/360) % 24, 1)\n", " \n", " if isinstance(datum, datetime):\n", " doys = int(datum.strftime('%j'))\n", " lts = calclt(datum, lon)\n", "\n", " elif isinstance(datum, list):\n", " doys = []\n", " lts = []\n", " for d, l in zip(datum, lon):\n", " doys.append(int(d.strftime('%j')))\n", " lts.append(calclt(d,l))\n", " \n", " return (doys, lts)\n", "\n", "\n", "def loadsatellitedata(starttime, endtime, satellite):\n", " '''Load data via viresclient\n", " \n", " Parameters\n", " ----------\n", " starttime : datetime\n", " endtime : datetime\n", " satellite : str\n", " name of satellite \n", " \n", " Returns\n", " -------\n", " pandas.DataFrame\n", " \n", " '''\n", " request = SwarmRequest()\n", " request.set_collection(f\"SW_OPER_IBI{satellite}TMS_2F\")\n", " request.set_products(measurements=request.available_measurements(\"IBI\"), \n", " auxiliaries=['F10_INDEX', 'F107'])\n", " data = request.get_between(starttime, endtime)\n", " df = data.as_dataframe()\n", " df = df[ (df.Bubble_Index > -1) ]\n", " \n", " return df\n", "\n", "\n", "def downloadfn(dates):\n", " '''Download Solar Radio Flux for specified date(s).\n", " \n", " Parameters\n", " ----------\n", " dates : datetime or list of datetimes\n", " \n", " Returns\n", " -------\n", " list of numpy.float64\n", " Solar Radio Flux\n", " '''\n", " def loaddf(year):\n", " df = pd.read_html(f'https://spaceweather.gc.ca/forecast-prevision/solar-solaire/solarflux/sx-5-flux-en.php?year={year}')[0]\n", " return df\n", " \n", " def getfn(df,date):\n", " return df[ (df.Date == date.strftime(\"%Y-%m-%d\")) & (df.Time == \"20:00:00\") ][\"Observed Flux\"].values[0]\n", " \n", " if isinstance(dates, datetime):\n", " years = [ dates.year ]\n", " dates = [ dates ]\n", " elif isinstance(dates, list):\n", " years = [ dates[0].year ]\n", " else: \n", " raise TypeError(f\"must be datetime or list of datetimes\")\n", " \n", " df_fn = loaddf(years[0])\n", " fn = []\n", " \n", " for d in dates:\n", " y = d.year\n", " if y not in years:\n", " df_fn = pd.concat([df_fn, loaddf(y)])\n", " fn.append(getfn(df_fn, d))\n", " years.append(y) \n", " else:\n", " fn.append(getfn(df_fn, d))\n", " \n", " return fn\n", "\n", "\n", "def calcindex(df):\n", " '''Calculation of the IBP index. If Solar Radio Flux lower than 80 the IBP index is set to -1.\n", " \n", " Parameters\n", " ----------\n", " df : pandas.DataFrame (with columns Doy, Lon, LT and F10.7)\n", " \n", " Returns\n", " -------\n", " pandas.DataFrame\n", " contains the columns: Doy (Day(s) of the year), Month (Month(s) from the day of the year), \n", " Lon (Longitude(s)), LT (Local Time(s)), F10.7 (solar index(es)), IBP (Ionospheric Bubble Index, value(s) between 0.0 and 1.0).\n", " '''\n", " def calcibp(index, row):\n", " df = ibp.calculateIBPindex(day_month=int(row[\"Doy\"]), longitude=row[\"Lon\"], local_time=row[\"LT\"], f107=row[\"F10.7\"])\n", " df = df.set_index([pd.Index([index])])\n", " return df\n", " \n", " index_sort = df.index\n", " \n", " tmp = df[ (df[\"F10.7\"] < 80) ].copy()\n", " tmp[\"IBP\"] = -1\n", " tmp['Month'] = [ int(datetime.strptime(f\"2023-{t}\",\"%Y-%j\").strftime('%m')) for t in tmp.Doy ]\n", " tmp['Month'] = tmp['Month'].astype(int)\n", " \n", " di = pd.DataFrame(columns=tmp.columns)\n", " \n", " i = 0\n", " for index, row in df[ (df[\"F10.7\"] >= 80) ].iterrows():\n", " if i == 0:\n", " di = calcibp(index, row)\n", " i = 1\n", " else:\n", " ib = calcibp(index, row)\n", " di = pd.concat([di, ib])\n", "\n", " di = pd.concat([di, tmp])\n", " di = di.reindex(index_sort)\n", " \n", " return di[['Doy','Month','Lon','LT','F10.7','IBP']]\n", "\n", "\n", "def satelliteIBP(df):\n", " '''Calculation of the IBP index at the specified time and longitude\n", " \n", " Parameters\n", " ----------\n", " df : pandas.DataFrame (with DatetimeIndex and columns Longitude and F107)\n", " \n", " Returns\n", " -------\n", " pandas.DataFrame\n", " with column IBP (Ionospheric Bubble Index, values between 0.0 and 1.0)\n", "\n", " '''\n", " \n", " # Create required columns\n", " df['Lon'] = np.rint(df['Longitude'])\n", " df['dt'] = df.index\n", " df['timediff'] = ((df.dt - df.dt.shift()) / np.timedelta64(1, 'm')).gt(3).cumsum()\n", " df['londiff'] = np.absolute(df['Lon'] - df['Longitude'])\n", "\n", " tmp = df.loc[df.groupby([\"timediff\", \"Lon\"]).agg({\"londiff\": ['idxmin']})[('londiff', 'idxmin')].values]\n", " tmp[\"Doy\"], tmp[\"LT\"] = getlt(tmp.index.tolist(), tmp.Lon)\n", " \n", " tmp.rename(columns={\"F107\": \"F10.7\"}, inplace=True)\n", " \n", " tmp['IBP'] = calcindex(tmp[['Doy', 'LT', 'Lon', 'F10.7']])['IBP']\n", " \n", " # merge df and tmp\n", " df = pd.concat([df, tmp['IBP']], axis=1)\n", " del df[\"Lon\"], df[\"dt\"], df[\"londiff\"], df[\"timediff\"]\n", " \n", " return df\n", "\n", "\n", "def globalIBP(datum, fn=150):\n", " '''Calculation of the IBP index for all longitudes at the specified time.\n", " \n", " Parameters\n", " ----------\n", " datum : datetime\n", " \n", " Returns\n", " -------\n", " pandas.DataFrame\n", " contains the columns: Doy (Day of the year), Month (Month from the day of the year), \n", " Lon (Longitudes), LT (Local Times), F10.7 (solar index), IBP (Ionospheric Bubble Index, values between 0.0 and 1.0).\n", " '''\n", " df = pd.DataFrame(columns=[\"Doy\", \"Lon\", \"LT\", \"F10.7\"])\n", " df[\"Lon\"] = np.arange(-180,180)\n", " df[\"Doy\"], df[\"LT\"] = getlt(datum, df[\"Lon\"])\n", " \n", " if datum.date() >= datetime.now().replace(tzinfo=timezone.utc).date():\n", " df[\"F10.7\"] = fn\n", " else:\n", " df[\"F10.7\"] = downloadfn(datum)[0]\n", " \n", " return calcindex(df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### IBP index around the earth at specific time" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "datum = datetime(2015, 1, 10, 2, 30).replace(tzinfo=timezone.utc)\n", "#datum = datetime(2016, 4, 21, 22,50,0).replace(tzinfo=timezone.utc)\n", "\n", "df1 = globalIBP(datum)\n", "df2 = globalIBP(datetime.now())\n", "\n", "fig, ax = plt.subplots(2,1,figsize=(12, 8), layout='constrained')\n", "df1.plot.scatter(x='Lon', y='IBP', title=f'IBP index at {datum.strftime(\"%Y-%m-%d %H:%M:%S\")}', \n", " xlabel='Longitude', ylabel='Probability', grid=True, xlim=(-180,179),\n", " c='IBP', colormap='plasma_r', vmin=0, vmax=1, ax=ax[0],\n", " sharex=False)\n", "# sharex=False - must be there because there is still a bug and xlabel and tickmarks will not appear otherwise\n", "df2.plot.scatter(x='Lon', y='IBP', title=f'IBP index now', \n", " xlabel='Longitude', ylabel='Probability', grid=True, xlim=(-180,179),\n", " c='IBP', colormap='plasma_r', vmin=0, vmax=1, ax=ax[1],\n", " sharex=False)\n", "\n", "for a in ax:\n", " a.set_ylim(top=1.05)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## IBP Index along the orbit of a satellite" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "starttime = datetime(2022, 11, 26, 17,0,0).replace(tzinfo=timezone.utc)\n", "endtime = datetime(2022, 11, 28, 5,0,0).replace(tzinfo=timezone.utc)\n", "\n", "sat = 'B'\n", "\n", "df = satelliteIBP(loadsatellitedata(starttime, endtime, sat))\n", "\n", "display(df)\n", "\n", "fig, axes = plt.subplots(1,2,figsize=(17, 5), layout='constrained')\n", "\n", "df.plot(x=\"Longitude\", y=\"Bubble_Index\", ax=axes[0], color='orange', kind='scatter', label='IBI - Bubble_Index')\n", "#df.plot(x=\"Longitude\", y=\"Bubble_Probability\", ax=axes[0], color='g', kind='scatter', label='IBI - Bubble_Probability')\n", "df.plot(x=\"Longitude\", y=\"IBP\", kind=\"scatter\", ax=axes[0], color='b', label='IBP')\n", "\n", "df.plot(y=\"Bubble_Index\", ax=axes[1], color='orange', marker=\".\", linewidth=0, label='IBI - Bubble_Index')\n", "#df.plot(x=\"Loy=\"Bubble_Probability\", ax=axes[0], color='g', kind='scatter', label='IBI - Bubble_Probability')\n", "df.plot(y=\"IBP\", marker=\".\", linewidth=0, ax=axes[1], color='b', label='IBP')\n", "\n", "\n", "axes[0].set_xlabel('Longitude')\n", "axes[0].set_xlim(-180,180)\n", "\n", "for ax in axes:\n", " ax.set_ylabel('probability')\n", " ax.set_ylim(top=1.05)\n", " ax.legend()\n", " ax.set_title(f'Plasma Bubbles detectet by satellite {sat}')\n", " ax.grid()\n", "\n", "plt.show() " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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" }, "vscode": { "interpreter": { "hash": "6a693262f012b0764ce1a509731aa7af644c1d221d30b43c6fed6082691c75d0" } } }, "nbformat": 4, "nbformat_minor": 2 }