satellite_elevation_azimuth_diagram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Satellite elevation-azimuth diagrams and PDOP\n",
    "### Guillaume Witz, Science IT Support, Bern University"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, re, glob\n",
    "from datetime import date, timedelta\n",
    "\n",
    "import numpy as np\n",
    "import scipy.signal\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import urllib.request\n",
    "from bz2 import BZ2File as bzopen\n",
    "\n",
    "import ipywidgets as ipw\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "\n",
    "import aiub"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#file location\n",
    "address = 'http://ftp.aiub.unibe.ch/BSWUSER52/STA/2019/'\n",
    "#address_sat = 'http://ftp.aiub.unibe.ch/CODE/'\n",
    "address_sat = 'http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/'\n",
    "#station coordinates file\n",
    "coord_file = 'COD19001.CRD.Z'\n",
    "#station to observe\n",
    "station_name = '14001M004'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elevation-azimuth\n",
    "#### Data import and computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import station data\n",
    "stations = aiub.import_stations(address,coord_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select station and calculate ellipsoidal coordinates\n",
    "curr_stat = stations[stations.statname == station_name].iloc[0]\n",
    "curr_stat['rad_stat_coord'] = aiub.cartesian_to_ellipsoidal(curr_stat.x_pos, curr_stat.y_pos,curr_stat.z_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2e727565d142423185ebcd3478ce33a6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(DatePicker(value=datetime.date(2018, 9, 1), description='Pick a start date', style=DescriptionS…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "date_min, date_max = aiub.date_picker()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import satellite data for all time-points\n",
    "temp_pd = aiub.import_sat_series(date_min.value, date_max.value, address_sat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#caluclate elevation and azimuth for all satellites for given station\n",
    "temp_pd = aiub.el_al_single_station_fast(temp_pd, curr_stat)\n",
    "temp_pd['curr_stat'] = curr_stat.statname"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd0d6873dfc84f8294998852bac4822b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value=\"<span style='float:left;font-size:2em;font-weight:bold'>Radial plot of satellite ob…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a3bdff26110a4495b16eea441cba8f55",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aiub.interactive_rad(temp_pd, stations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PDOP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### User input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bf201f18f06843e99402946fb0753951",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(DatePicker(value=datetime.date(2018, 9, 1), description='Pick a start date', style=DescriptionS…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "date_min, date_max = aiub.date_picker()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data import and calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a list of dates between start and end date \n",
    "date_list = [date_min.value + timedelta(days=x) for x in range((date_max.value-date_min.value).days+1)]\n",
    "#calculate gps weeks times (year, week, day)\n",
    "gps_weeks = [aiub.date_to_gpsweeks(x) for x in date_list]\n",
    "\n",
    "#load satellite data and calculate el, az\n",
    "temp_pd = pd.concat([aiub.import_RM_file(address_sat, g) for g in gps_weeks]).reset_index()\n",
    "temp_pd = aiub.el_al_single_station_fast(temp_pd, curr_stat)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove elevation < 5°\n",
    "temp_pd = temp_pd[temp_pd.e >2*np.pi*5/360]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate norm of satellite-station vector Xgeo\n",
    "temp_pd['Xgeo_norm'] = temp_pd.apply(lambda row: np.linalg.norm([row['Xgeo'],row['Ygeo'],row['Zgeo']]),axis = 1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_grouped = temp_pd.set_index('datetime').groupby(pd.Grouper(freq='30min'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/gw18g940/miniconda3/envs/astronomy/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in sqrt\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "#calculate pdop for each group\n",
    "pdop = dict((key, []) for key in temp_pd.satellite.unique())\n",
    "\n",
    "for name, df0 in time_grouped:\n",
    "    sat_grouped = df0.groupby('satellite')\n",
    "    for name2, df in sat_grouped:\n",
    "        A_mat = np.stack([-df['Xgeo']/df['Xgeo_norm'],-df['Ygeo']/df['Xgeo_norm'],\n",
    "              -df['Zgeo']/df['Xgeo_norm'],np.ones(len(df))],axis = 1)\n",
    "\n",
    "        try:\n",
    "            inv_mat = np.linalg.inv(np.matmul(A_mat.T,A_mat))\n",
    "        except:\n",
    "            False\n",
    "\n",
    "        pdop[name2].append(np.sqrt(inv_mat[0,0]+inv_mat[1,1]+inv_mat[2,2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate power spectrum and set the right (?) scale \n",
    "powerspec = {}\n",
    "for x in pdop.keys():\n",
    "    powerspec[x] = scipy.signal.periodogram(pdop[x],fs = 48, scaling='spectrum')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define interactive plotting function \n",
    "def plot_pdop(sat_type, pdop, powerspec):\n",
    "    fig, ax = plt.subplots(2,1)\n",
    "    for s in sat_type:\n",
    "        ax[0].plot(np.arange(len(pdop[s]))/(24*2),pdop[s], label = s)\n",
    "        \n",
    "        ax[1].plot(powerspec[s][0],powerspec[s][1],'-o', label = s)\n",
    "    ax[0].legend()\n",
    "    \n",
    "    #ax[1].set_xlim\n",
    "    plt.xscale('log')\n",
    "    ax[1].set_xlim((0,10))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5e7dbc270e8542b6995587fcb5359e86",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(SelectMultiple(description='sat_type', index=(0,), options=('G', 'R', 'E', 'C'), value=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ipw.interactive(plot_pdop, sat_type = ipw.SelectMultiple(options =temp_pd.satellite.unique(),index=[0]),\n",
    "                pdop = ipw.fixed(pdop), powerspec = ipw.fixed(powerspec))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}