ant_pos_altaz.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ccff02-7bbf-4e09-9b35-f7eb44ef892c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3b5be7f8-2a8b-416d-9d63-5097cc34adee",
   "metadata": {},
   "source": [
    "%%time\n",
    "\n",
    "from pathlib import Path\n",
    "from collections import namedtuple\n",
    "\n",
    "from astropy.table import Table\n",
    "from astropy.time import Time\n",
    "from astropy.coordinates import Angle, SkyCoord, EarthLocation, AltAz\n",
    "import astropy.units as u\n",
    "import geoviews\n",
    "import pytz\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import datashader\n",
    "\n",
    "import holoviews\n",
    "import geopandas\n",
    "import geoviews\n",
    "from cartopy import crs\n",
    "import colorcet as cc\n",
    "import geoviews\n",
    "import seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "geoviews.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f47a5d4a-de34-4191-8db9-e97753bdb0e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# filename = \"AGBT16A_353.AltAz.parquet\"\n",
    "filename = \"ant_pos_all_v2.parquet\"\n",
    "# filename = \"ant_pos_with_az_al_small.parquet\"\n",
    "df = pd.read_parquet(filename, columns=[\"DMJD\", \"RAJ2000\", \"DECJ2000\", \"MNT_AZ\", \"MNT_EL\"])\n",
    "lat_col = \"MNT_EL\"\n",
    "lon_col = \"MNT_AZ\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd17a26-5cfc-4935-9d43-df664b68a968",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scs = SkyCoord(df[\"RAJ2000\"], df[\"DECJ2000\"], equinox=\"J2000\", unit=u.degree)\n",
    "\n",
    "# gbt = EarthLocation(lat=38.433119 * u.deg, lon=79.839833 * u.deg, height=824.595 * u.m)\n",
    "# frame = AltAz(obstime=Time(df[\"DMJD\"], format=\"mjd\"), location=gbt)\n",
    "# altaz = scs.transform_to(frame).to_table()\n",
    "# df[\"beam_az\"] = altaz[\"az\"]\n",
    "# df[\"beam_el\"] = altaz[\"alt\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bff6b44-633e-45b0-bde7-1276ac2dab9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True)\n",
    "seaborn.boxplot(x=df[\"MNT_AZ\"], ax=ax1)\n",
    "seaborn.boxplot(x=df[\"MNT_EL\"], ax=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2fbe334-8f72-4ba6-995e-9fd23384903b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# wrap\n",
    "df[df[\"MNT_AZ\"] >= 360] = df[df[\"MNT_AZ\"] >= 360] - 360\n",
    "df[df[\"MNT_AZ\"] < 0] = df[df[\"MNT_AZ\"] < 0] + 360"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45aeb5be-0b4c-438e-90a9-b1f58a713b03",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter out outliers\n",
    "df = df[(df[\"MNT_AZ\"] < 360) & (df[\"MNT_EL\"] < 90)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba86e18e-173d-43d6-8fb9-2d60da054548",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[(df[\"MNT_EL\"] > 0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc347403-7063-4025-81f9-14aa5af2870b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True)\n",
    "seaborn.boxplot(x=df[\"MNT_AZ\"], ax=ax1)\n",
    "seaborn.boxplot(x=df[\"MNT_EL\"], ax=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d342bc1b-cb6d-47a5-90b3-9588dc9afe84",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "points = geoviews.Points(df, [lon_col, lat_col])\n",
    "projection = crs.AzimuthalEquidistant(central_longitude=180, central_latitude=90, false_northing=0)\n",
    "# projection = crs.AzimuthalEquidistant(central_longitude=79.839833, central_latitude=90-38.433119, false_northing=0)\n",
    "# projection = crs.AzimuthalEquidistant(central_longitude=0, central_latitude=38.433119, false_northing=0)\n",
    "projected = geoviews.operation.project_points(points, projection=projection)\n",
    "start = Time(projected.data[projected.data[\"DMJD\"] > 50_000][\"DMJD\"].min(), format=\"mjd\").to_datetime(timezone=pytz.UTC)\n",
    "end = Time(projected.data[\"DMJD\"].max(), format=\"mjd\").to_datetime(timezone=pytz.UTC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a49673-7146-40d3-9e8d-ab39a4e35581",
   "metadata": {},
   "outputs": [],
   "source": [
    "label_locs = pd.DataFrame([(v, 15 if v >= 0 else -15, str(v)) for v in range(-180, 180, 30)], columns=[\"lat\", \"lon\", \"text\"])\n",
    "print(label_locs)\n",
    "projected_labels = geoviews.operation.project_points(geoviews.Points(label_locs, [\"lon\", \"lat\"]), projection=projection)\n",
    "\n",
    "projected_labels.data\n",
    "\n",
    "labels = holoviews.Labels({('x', 'y'): projected_labels.data[[\"lon\", \"lat\"]], 'text': projected_labels.data[\"text\"]}, ['x', 'y'], 'text')\n",
    "# projected_labels[[\"lon\", \"lat\"]]\n",
    "print(projected_labels.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f84dfbd2-a89c-42f0-9bbf-1b8e8981d0b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "def plot_az_el(projected, width=800):\n",
    "    hv_shaded = holoviews.operation.datashader.datashade(projected, dynamic=True, precompute=True, cmap=cc.rainbow, x_range=(-10e6, 10e6), y_range=(-10e6, 10e6))\n",
    "    hv_spread = holoviews.operation.datashader.dynspread(hv_shaded, threshold=0.1, max_px=1, dynamic=True)\n",
    "    hv_spread = hv_spread.opts(bgcolor=\"black\", width=width, height=width, global_extent=False, projection=projection)\n",
    "    return hv_spread * geoviews.feature.grid\n",
    "    \n",
    "\n",
    "from datetime import datetime\n",
    "from tqdm import tqdm\n",
    "geoviews.extension('bokeh')\n",
    "plots = []\n",
    "total_width=4000\n",
    "num_cols=5\n",
    "# Skip 2002 since it is only partially complete\n",
    "for year in tqdm(range(start.year + 1, end.year + 1)):\n",
    "    start_mjd = Time(datetime(year, 1, 1)).to_value(\"mjd\")\n",
    "    end_mjd = Time(datetime(year + 1, 1, 1)).to_value(\"mjd\")\n",
    "    data_for_year = projected.data[(projected.data[\"DMJD\"] >= start_mjd) & (projected.data[\"DMJD\"] < end_mjd)]\n",
    "    \n",
    "    plots.append(plot_az_el(geoviews.Points(data_for_year, crs=projection), width=total_width//num_cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5206a4d2-de5e-4db7-b147-6689f76b2137",
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = holoviews.core.Layout(plots).cols(num_cols)\n",
    "layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3750565-272b-4abc-8c2c-8b17c129c2bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# output_x = 2000\n",
    "# hv_shaded = holoviews.operation.datashader.datashade(projected, dynamic=False, precompute=True, cmap=cc.rainbow, x_range=(-10e6, 10e6), y_range=(-10e6, 10e6))\n",
    "# hv_spread = holoviews.operation.datashader.dynspread(hv_shaded, threshold=0.4, max_px=1)\n",
    "# hv_spread = hv_spread.opts(bgcolor=\"white\", width=output_x, height=output_x, global_extent=False, projection=projection)\n",
    "# grid = geoviews.feature.grid\n",
    "# # hv_spread * grid * labels\n",
    "holoviews.save(layout, \"altaz.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b3b7c7-a0c6-4115-9fb8-aa58be13926c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# geoviews.extension('matplotlib')\n",
    "\n",
    "# from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "# from datashader.mpl_ext import dsshow, alpha_colormap\n",
    "\n",
    "# fig = plt.figure(figsize=(9, 9))\n",
    "\n",
    "# grid = ImageGrid(fig, 111, nrows_ncols=(2, 2), axes_pad=0.5, share_all=True, cbar_location=\"right\", cbar_mode=\"each\", cbar_size=\"5%\", cbar_pad=\"2%\")\n",
    "\n",
    "# artist0 = dsshow(df, datashader.Point('MNT_AZ', 'MNT_EL'), datashader.count(), vmax=1000, aspect='equal', ax=grid[0])"
   ]
  }
 ],
 "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",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}