dva_dwingeloo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "649a2e98",
   "metadata": {},
   "source": [
    "### Import necessary packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9cf2132",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from astropy.wcs import WCS\n",
    "from reproject import reproject_interp\n",
    "import numpy as np\n",
    "import astropy.io.fits as fits\n",
    "import healpy as hp\n",
    "from matplotlib import cm\n",
    "from astropy_healpix import HEALPix\n",
    "from astropy.coordinates import SkyCoord, ICRS, Galactic\n",
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4275f57",
   "metadata": {},
   "source": [
    "### Read in Dwingeloo data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aedd8bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "hduQ = fits.open('/srv/data/dwingeloo/DWING_820_EQU_Q.HPX.fits')\n",
    "hdrQ = hduQ[0].header\n",
    "dataQ = hduQ[0].data[0]\n",
    "print(dataQ.shape)\n",
    "\n",
    "hduU = fits.open('/srv/data/dwingeloo/DWING_820_EQU_U.HPX.fits')\n",
    "hdrU = hduU[0].header\n",
    "dataU = hduU[0].data[0]\n",
    "print(dataU.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da41e03d",
   "metadata": {},
   "source": [
    "### Pick the range of RA and dec to plot, and the gridding pixel size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f7ec79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample a grid in RA/Dec\n",
    "RA_min = 0\n",
    "RA_max = 360\n",
    "dec_min = -10\n",
    "dec_max = 90\n",
    "delta = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43d70b4c",
   "metadata": {},
   "source": [
    "### Set up the HEALpix projection, grid and coordinates\n",
    "#### NOTE: if all files have the same NSIDE, just need to do this once"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c60e72c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Healpix parameters:\n",
    "nside = hdrQ['NSIDE']\n",
    "order = hdrQ['ORDERING']\n",
    "hpx = HEALPix(nside=nside, order=order, frame=ICRS())\n",
    "\n",
    "# Number of pixels along each axis:\n",
    "numRA = int((RA_max-RA_min)/delta)\n",
    "numdec = int((dec_max-dec_min)/delta)\n",
    "\n",
    "# 1D and 2D RA and dec grids:\n",
    "ra = np.linspace(RA_min, RA_max, numRA) * u.deg\n",
    "dec = np.linspace(dec_min, dec_max, numdec) * u.deg\n",
    "ra_grid, dec_grid = np.meshgrid(ra, dec)\n",
    "\n",
    "# Set up Astropy coordinate objects\n",
    "coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa7dfd44",
   "metadata": {},
   "source": [
    "### Reproject the HEALpix values onto the 2D grid for each map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270be18b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Qmap = hpx.interpolate_bilinear_skycoord(coords, dataQ)\n",
    "Qmap = Qmap.reshape((numdec, numRA))[:,::-1]\n",
    "print(Qmap.shape)\n",
    "\n",
    "Umap = hpx.interpolate_bilinear_skycoord(coords, dataU)\n",
    "Umap = Umap.reshape((numdec, numRA))[:,::-1]\n",
    "print(Umap.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8644f1e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdr_new = fits.Header.fromstring(\"\"\"\n",
    "NAXIS   =                    2\n",
    "CUNIT1  = 'deg     '\n",
    "CUNIT2  = 'deg     '\n",
    "\"\"\", sep='\\n')\n",
    "\n",
    "hdr_new['COORDSYS']  = 'icrs'\n",
    "hdr_new['NAXIS1']  = Qmap.shape[1] \n",
    "hdr_new['NAXIS2']  = Qmap.shape[0]\n",
    "\n",
    "hdr_new['CTYPE1']  = 'RA---CAR'\n",
    "hdr_new['CRPIX1']  = Qmap.shape[1]/2.+0.5 \n",
    "hdr_new['CRVAL1']  = (RA_min+RA_max)/2.          \n",
    "hdr_new['CDELT1']  = -delta\n",
    "\n",
    "hdr_new['CTYPE2']  = 'DEC--CAR'\n",
    "hdr_new['CRPIX2']  = Qmap.shape[0]/2.+0.5\n",
    "hdr_new['CRVAL2']  = (dec_min+dec_max)/2.\n",
    "hdr_new['CDELT2']  = delta\n",
    "\n",
    "print(repr(hdr_new))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6578e72a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fits.writeto('/home/aordog/DATA/fits_map_test.fits',Qmap,hdr_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e2f583e",
   "metadata": {},
   "source": [
    "### Plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37bbedfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 4))\n",
    "im = plt.imshow(Qmap, extent=[RA_max, RA_min, dec_min, dec_max], cmap=plt.cm.RdBu_r, \n",
    "                origin='lower', aspect='equal',vmin=-1,vmax=1)\n",
    "plt.colorbar(im)\n",
    "plt.xlabel('Right ascension (deg)')\n",
    "plt.ylabel('Declination (deg)')\n",
    "plt.title('Stokes Q')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(20, 4))\n",
    "im = plt.imshow(Umap, extent=[RA_max, RA_min, dec_min, dec_max], cmap=plt.cm.RdBu_r, \n",
    "                origin='lower', aspect='equal',vmin=-1,vmax=1)\n",
    "plt.colorbar(im)\n",
    "plt.xlabel('Right ascension (deg)')\n",
    "plt.ylabel('Declination (deg)')\n",
    "plt.title('Stokes U')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfcd4b14",
   "metadata": {},
   "outputs": [],
   "source": [
    "l_list = []\n",
    "b_list = []\n",
    "ra_list = []\n",
    "dec_list = []\n",
    "PI_list = []\n",
    "gal_theta_list = []\n",
    "eq_theta_list = []\n",
    "with open('/srv/data/dwingeloo_data/original/out820') as f:\n",
    "    for line in f:\n",
    "        line = f.readline()\n",
    "        l_list.append(float(line[0:7]))\n",
    "        b_list.append(float(line[7:13]))\n",
    "        ra_list.append(float(line[13:19]))\n",
    "        dec_list.append(float(line[19:25]))\n",
    "        PI_list.append(float(line[25:31]))\n",
    "        gal_theta_list.append(float(line[31:37]))\n",
    "        eq_theta_list.append(float(line[37:43]))\n",
    "f.close()\n",
    "ra_arr = np.asarray(ra_list)\n",
    "dec_arr = np.asarray(dec_list)\n",
    "\n",
    "dec_arr2 = dec_arr.copy()\n",
    "dec_arr2[np.where(dec_arr>=270.)] = dec_arr2[np.where(dec_arr>=270.)]-360.\n",
    "\n",
    "Qpts = np.asarray(PI_list)*np.cos((2*np.asarray(eq_theta_list))*np.pi/180)\n",
    "Upts = np.asarray(PI_list)*np.sin((2*np.asarray(eq_theta_list))*np.pi/180)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6048bbe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 6))\n",
    "im = plt.imshow(Qmap, extent=[RA_max, RA_min, dec_min, dec_max], cmap=plt.cm.RdBu_r, \n",
    "                origin='lower', aspect='equal',vmin=-1,vmax=1)\n",
    "plt.scatter(ra_arr,dec_arr2,c=Qpts,vmin=-1,vmax=1,cmap=plt.cm.RdBu_r,s=30,edgecolors='k')\n",
    "plt.colorbar(im)\n",
    "plt.xlabel('Right ascension (deg)')\n",
    "plt.ylabel('Declination (deg)')\n",
    "plt.title('Stokes Q')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(20, 6))\n",
    "im = plt.imshow(Umap, extent=[RA_max, RA_min, dec_min, dec_max], cmap=plt.cm.RdBu_r, \n",
    "                origin='lower', aspect='equal',vmin=-1,vmax=1)\n",
    "plt.scatter(ra_arr,dec_arr2,c=Upts,vmin=-1,vmax=1,cmap=plt.cm.RdBu_r,s=30,edgecolors='k')\n",
    "plt.colorbar(im)\n",
    "plt.xlabel('Right ascension (deg)')\n",
    "plt.ylabel('Declination (deg)')\n",
    "plt.title('Stokes U')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74b0fbf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 6))\n",
    "im = plt.imshow(0.5*np.arctan2(Umap,Qmap)*180./np.pi, extent=[RA_max, RA_min, dec_min, dec_max], cmap=plt.cm.hsv, \n",
    "                origin='lower', aspect='equal',vmin=-90,vmax=90)\n",
    "#plt.scatter(ra_arr,dec_arr2,c=0.5*np.arctan2(Upts,Qpts)*180./np.pi,vmin=-90,vmax=90,cmap=plt.cm.hsv,s=40,edgecolors='k')\n",
    "plt.colorbar(im)\n",
    "plt.xlabel('Right ascension (deg)')\n",
    "plt.ylabel('Declination (deg)')\n",
    "plt.title('PA')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3b39656",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(eq_theta_list,0.5*np.arctan2(Upts,Qpts)*180./np.pi)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd2e8177",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import angular_separation\n",
    "\n",
    "qwcs = WCS(hdr_new)\n",
    "\n",
    "print(Qmap.shape[1])\n",
    "print(Qmap.shape[0])\n",
    "print(repr(qwcs))\n",
    "\n",
    "c = SkyCoord.from_pixel(np.arange(Qmap.shape[1])[:,np.newaxis], np.arange(Qmap.shape[0])[np.newaxis,:], qwcs.celestial)\n",
    "radius = 1 #* u.deg\n",
    "drao_q_avg = []\n",
    "drao_u_avg = []\n",
    "for i in np.arange(len(ra_arr)):\n",
    "    #print(ra_arr[i])\n",
    "    #print(c.icrs.ra.deg)\n",
    "    sep = angular_separation(np.radians(ra_arr[i]), np.radians(dec_arr2[i]), c.icrs.ra.radian, c.icrs.dec.radian)*180/np.pi\n",
    "    #print(type(sep))\n",
    "    idx = np.where(sep <= radius)\n",
    "    #print(i)\n",
    "    #print(ra_arr[i],dec_arr2[i])\n",
    "    #print(c.icrs.ra.deg[idx])\n",
    "    #print(c.icrs.dec.deg[idx])\n",
    "    #print('')\n",
    "    drao_q_avg.append(np.nanmean(Qmap[idx[1], idx[0]]))\n",
    "    drao_u_avg.append(np.nanmean(Umap[idx[1], idx[0]]))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d9e51a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(drao_q_avg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "822733ef",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}