airmass.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "%pylab inline\n",
    "from scipy import stats\n",
    "from angles import r2d, r2arcs, d2arcs, arcs2r\n",
    "def arcm2r(theta):\n",
    "    return arcs2r(theta*60)\n",
    "import seaborn as sns;sns.set_style('darkgrid')\n",
    "import lsst.sims.maf.stackers as stackers\n",
    "import treecorr as tr\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "matplotlib.rcParams['figure.figsize'] = 19,14\n",
    "import lsst.sims.maf.slicers as slicers\n",
    "import lsst.sims.maf.metrics as metrics\n",
    "import lsst.sims.maf.metricBundles as metricBundles\n",
    "import lsst.sims.maf.db as db\n",
    "from collections import defaultdict\n",
    "import treecorr\n",
    "from itertools import chain\n",
    "sns.set_context('poster',font_scale=1.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outDir='temp'\n",
    "conn = db.OpsimDatabaseV4('simruns/feature_baseline_update_10yrsv2')\n",
    "bundleList = []\n",
    "sql = 'filter = \"r\" and proposalId = 0 and airmass>0'\n",
    "metric = metrics.CountMetric('airmass')\n",
    "slicer = slicers.HealpixSlicer()\n",
    "#dith_stack = stackers.NormAirmassStacker()\n",
    "bundleList.append(metricBundles.MetricBundle(metric, slicer, sql))#, stackerList=[dith_stack]))\n",
    "bd = metricBundles.makeBundlesDictFromList(bundleList)\n",
    "mbg = metricBundles.MetricBundleGroup(bd, conn, outDir=outDir, saveEarly= False)\n",
    "mbg.runAll()\n",
    "conn.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num, xedges, yedges = np.histogram2d(mbg.simData['fieldRA'], mbg.simData['fieldDec'], bins=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H, xedges, yedges = np.histogram2d(mbg.simData['fieldRA'], mbg.simData['fieldDec'], weights=mbg.simData['airmass'], bins=15)\n",
    "H = H.T\n",
    "plt.imshow(H/num.T, interpolation='cubic', origin='low', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- query the database for filter, dec, airmass\n",
    "- for each filter: plot distribution of airmasses histogram overall\n",
    "- for each filter: for each declination band: plot airmass histogram\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outDir='temp'\n",
    "conn = db.OpsimDatabaseV4('simruns/feature_baseline_update_10yrsv2')\n",
    "bundleList = []\n",
    "sql = 'proposalId = 0 and airmass>0'\n",
    "metric = metrics.CountMetric('filter')\n",
    "slicer = slicers.HealpixSlicer()\n",
    "dith_stack = stackers.NormAirmassStacker()\n",
    "bundleList.append(metricBundles.MetricBundle(metric, slicer, sql, stackerList=[dith_stack]))\n",
    "bd = metricBundles.makeBundlesDictFromList(bundleList)\n",
    "mbg = metricBundles.MetricBundleGroup(bd, conn, outDir=outDir, saveEarly= False)\n",
    "mbg.runAll()\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_data = mbg.simData['filter']\n",
    "dec_data = mbg.simData['fieldDec']\n",
    "airmass_data = mbg.simData['airmass']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "ks_tuple = namedtuple('ks_tuple', ['Dstatistic', 'pvalue'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airmass_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks_stats_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_airmass = airmass_data/min(airmass_data)\n",
    "filters = ['u','g','r','i','z','y']\n",
    "ks_stats_dict = {k:[] for k in filters}\n",
    "airmass_dict = {k:[] for k in filters}\n",
    "for dec_band, dec_band_next in zip([0,-20,-40,-60,-80],[20,0,-20,-40,-60]):\n",
    "    airmass_per_dec_band = airmass_data[np.logical_and(dec_data>dec_band, dec_data<dec_band_next)]\n",
    "    plt.hist(airmass_per_dec_band, label='average', alpha=0.2, bins=100)\n",
    "    for filt in filters:\n",
    "        \n",
    "        plt.axvline(2)\n",
    "        plt.axvline(1)\n",
    "        filtered_airmass = airmass_data[np.logical_and(np.logical_and(filter_data==filt,dec_data>dec_band),dec_data<dec_band_next)]\n",
    "        #filtered_airmass /= min(filtered_airmass)\n",
    "        plt.hist(filtered_airmass, label=filt+str(dec_band), histtype='step', linewidth=2.5, bins=100)\n",
    "        airmass_dict[filt].append(filtered_airmass)\n",
    "        ks_test = stats.ks_2samp(filtered_airmass, airmass_per_dec_band)\n",
    "        ks_stats_dict[filt].append(ks_tuple(Dstatistic=ks_test.statistic, pvalue=ks_test.pvalue))\n",
    "    plt.legend()\n",
    "    plt.xscale('log')\n",
    "    plt.xlim(0,2)\n",
    "    plt.yscale('log')\n",
    "    plt.figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "ra_data = mbg.simData['fieldRA']\n",
    "filters = ['u','g','r','i','z','y']\n",
    "dec_bands = [0,-10,-20,-30,-40,-50,-60,-70]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks_stats_dict = {k:{dec:[] for dec in dec_bands} for k in filters}\n",
    "for dec_band, dec_band_next in zip([0,-10,-20,-30,-40,-50,-60,-70],\n",
    "                                   [10,0,-10,-20,-30,-40,-50,-60]):\n",
    "    for filt in filters:\n",
    "        airmass_perfilter_perdec = airmass_data[reduce(np.logical_and,\n",
    "                                                       (filter_data==filt,\n",
    "                                                        dec_data>dec_band,\n",
    "                                                        dec_data<dec_band_next))]\n",
    "        for ra_band, ra_band_next in zip(np.linspace(0,350,36),np.linspace(10,360,36)):\n",
    "            filtered_airmass = airmass_data[reduce(np.logical_and,\n",
    "                                                   (filter_data==filt,\n",
    "                                                    dec_data>dec_band,\n",
    "                                                    dec_data<dec_band_next, \n",
    "                                                    ra_data>ra_band, \n",
    "                                                    ra_data<ra_band_next))]\n",
    "            ks_test = stats.ks_2samp(filtered_airmass, airmass_perfilter_perdec)\n",
    "            if ks_test.pvalue<10**-10:\n",
    "                pval = 10**-10\n",
    "            else:\n",
    "                pval = ks_test.pvalue\n",
    "            ks_stats_dict[filt][dec_band].append(ks_tuple(Dstatistic=ks_test.statistic, \n",
    "                                                          pvalue=pval))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks_stats_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dec in dec_bands:\n",
    "    plt.figure()\n",
    "    plt.title(str(dec))\n",
    "    for filt in filters:\n",
    "        pvals = [i.pvalue for i in ks_stats_dict[filt][dec]]\n",
    "        plt.hist(pvals, alpha=0.2, linewidth=2.5, label=filt, bins=np.logspace(-10,0,20))\n",
    "        plt.xscale('log')\n",
    "        plt.legend()\n",
    "        plt.gca().set_xlim(right=2)\n",
    "        plt.yscale('log')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lsst",
   "language": "python",
   "name": "lsst"
  },
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}