star_density.py

import numpy as np
from astroquery.vizier import Vizier
import pylab as pl
from astropy import constants, units as u, table, stats, coordinates, wcs, log, coordinates as coord, convolution, modeling, time; from astropy.io import fits, ascii

Vizier.ROW_LIMIT = 3e5
rslt = Vizier.query_constraints(catalog="II/246", GLAT='<5 & >-5', Kmag='<7', GLON='<45 | >315')[0]
crds = coordinates.SkyCoord(rslt['RAJ2000'], rslt['DEJ2000'], frame='fk5', unit=(u.deg, u.deg)).galactic

gridy,gridx = np.mgrid[-5:5:0.333, -45:45:0.333]
his,xe,ye = np.histogram2d(crds.galactic.l.wrap_at(180*u.deg).deg, crds.galactic.b.deg, bins=[gridx[0,:], gridy[:,0]])

fig = pl.figure(1)
fig.clf()
ax = fig.gca()
pl.imshow(his.T, extent=[xe.min(),xe.max(),ye.min(),ye.max()])
ax.set_aspect('equal') #gridx.shape[1]/gridx.shape[0])
cax = fig.add_axes([0.93, 0.4, 0.02, 0.2])
pl.colorbar(cax=cax)

ax.set_xlabel("Galactic Longitude")
ax.set_ylabel("Galactic Latitude")
ax.set_title("Number of sources with K<7 per 20x20 arcminute region")
pl.savefig("galactic_klt7_density.png", bbox_inches='tight')

for kmax in range(3,7):
    kltX = rslt['Kmag'] < kmax

    gridy,gridx = np.mgrid[-5:5:0.333, -45:45:0.333]
    his,xe,ye = np.histogram2d(crds[kltX].galactic.l.wrap_at(180*u.deg).deg,
                               crds[kltX].galactic.b.deg,
                               bins=[gridx[0,:], gridy[:,0]])

    fig.clf()
    ax = fig.gca()
    pl.imshow(his.T, extent=[xe.min(),xe.max(),ye.min(),ye.max()])
    ax.set_aspect('equal') #gridx.shape[1]/gridx.shape[0])
    cax = fig.add_axes([0.93, 0.4, 0.02, 0.2])
    pl.colorbar(cax=cax)

    ax.set_xlabel("Galactic Longitude")
    ax.set_ylabel("Galactic Latitude")
    ax.set_title(f"Number of sources with K<{kmax} per 20x20 arcminute region")

    pl.savefig(f"galactic_klt{kmax}_density.png", bbox_inches='tight')