########################################################
# Started Logging At: 2024-10-08 08:46:03
########################################################
########################################################
# # Started Logging At: 2024-10-08 08:46:03
########################################################
from astropy import table
sanders = Table.read('../sgrE.csv')
from astropy import table
from astropy.table import Table
sanders = Table.read('../sgrE.csv')
sanders
#[Out]# <Table length=197671>
#[Out]#    sourceid    astfit_epochs ...        epmb             pml_pmb_corr    
#[Out]#     int64          int64     ...      float64              float64       
#[Out]# -------------- ------------- ... ------------------ ---------------------
#[Out]# 13640652012519           131 ...  3.796065772319169  0.040963077269144645
#[Out]# 13640652021185           199 ...  0.823237672712744   0.05732593534318079
#[Out]# 13640652002403            36 ...  9.439942894239984   0.02061710356674774
#[Out]# 13640652015378           201 ... 0.3160578479498499   0.03417549108494654
#[Out]# 13644749012214           168 ... 1.2819993988827965   0.07018086251265186
#[Out]# 13644749022409           180 ... 1.2734729601558754   0.05104081017850679
#[Out]# 13644749011280           128 ...  4.037744856304961  0.043107124049581697
#[Out]# 13644749014919           185 ... 1.4272258491882222   0.04790572495368927
#[Out]# 13644749012899            70 ...  5.030027447479226  0.041335465334242615
#[Out]#            ...           ... ...                ...                   ...
#[Out]# 13611984022351            15 ...  7.227753076454498  -0.10370006465266776
#[Out]# 13611984011380           239 ... 0.4818779058897811   -0.1442519256428744
#[Out]# 13611984006606           221 ...   2.18524965088116  -0.09655196661681018
#[Out]# 13611984009283            72 ...  4.107072247171149  -0.10576486770587286
#[Out]# 13611984013410           205 ... 1.7150833638611107  -0.09258263053844404
#[Out]# 13611984006630           208 ... 2.5090553997059954  -0.09388323269336052
#[Out]# 13611984005156           228 ... 0.8599608461805929  -0.03609533824754119
#[Out]# 13611984010955           258 ... 0.7395392251943154 -0.021174377312691093
#[Out]# 13611984016555            12 ... 12.295529517934133     -0.29187669234933
#[Out]# 13616080021773           125 ...  5.092177456309112   -0.1250389363841114
import pylab as pl
from astropy import table
from astropy.table import Table
sanders = Table.read('../sgrE.csv')
significant = (sanders['pmra'] > sanders['pmra_error'] * 2) & (sanders['pmdec'] > sanders['pmdec_error'] * 2)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], size=1, alpha=0.5)
significant = (sanders['pmra'] > sanders['pmra_error'] * 2) & (sanders['pmdec'] > sanders['pmdec_error'] * 2)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
#[Out]# <matplotlib.collections.PathCollection at 0x152ac6c85190>
significant = ((sanders['pmra'] > sanders['pmra_error'] * 2) &
               (sanders['pmdec'] > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.
significant = ((sanders['pmra'] > sanders['pmra_error'] * 2) &
               (sanders['pmdec'] > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
#[Out]# <matplotlib.collections.PathCollection at 0x152ac64fc9b0>
significant = ((sanders['pmra'] > sanders['pmra_error'] * 2) &
               (sanders['pmdec'] > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
significant = ((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
               (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
import pylab as pl
import numpy as np
from astropy import table
from astropy.table import Table
sanders = Table.read('../sgrE.csv')
sanders[:3]
#[Out]# <Table length=3>
#[Out]#    sourceid    astfit_epochs ...        epmb           pml_pmb_corr    
#[Out]#     int64          int64     ...      float64            float64       
#[Out]# -------------- ------------- ... ----------------- --------------------
#[Out]# 13640652012519           131 ... 3.796065772319169 0.040963077269144645
#[Out]# 13640652021185           199 ... 0.823237672712744  0.05732593534318079
#[Out]# 13640652002403            36 ... 9.439942894239984  0.02061710356674774
significant = ((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
               (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
significant = ((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
               (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
               (sanders['pmra_error'] < 10) &
               (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.5)
pl.axis([-30, 30, -30, 30]);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.1)
pl.axis([-30, 30, -30, 30]);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.axis([-30, 30, -30, 30]);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.axis([-15, 15, -15, 15]);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_k_mean_mag'], sanders['phot_k_mean_mag'])
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'], sanders['phot_ks_mean_mag'])
#[Out]# (array([[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
#[Out]#          0.0000e+00, 0.0000e+00, 3.0000e+00, 3.0000e+00, 0.0000e+00],
#[Out]#         [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
#[Out]#          0.0000e+00, 0.0000e+00, 1.0000e+00, 6.0000e+00, 1.0000e+00],
#[Out]#         [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
#[Out]#          0.0000e+00, 3.0000e+00, 1.2000e+01, 4.7000e+01, 1.0000e+01],
#[Out]#         [1.0000e+00, 2.6000e+01, 1.0300e+02, 2.1100e+02, 5.6500e+02,
#[Out]#          1.1400e+03, 2.0740e+03, 4.8950e+03, 1.9117e+04, 1.2638e+04],
#[Out]#         [2.0000e+00, 1.3000e+02, 6.4100e+02, 1.3930e+03, 2.6030e+03,
#[Out]#          8.4500e+03, 8.4170e+03, 1.2301e+04, 2.2391e+04, 4.2970e+03],
#[Out]#         [1.6000e+01, 6.7600e+02, 1.7040e+03, 2.9610e+03, 5.7510e+03,
#[Out]#          1.4746e+04, 2.6263e+04, 1.8354e+04, 1.0433e+04, 2.2400e+02],
#[Out]#         [2.0000e+00, 1.5900e+02, 4.1900e+02, 6.2100e+02, 9.6100e+02,
#[Out]#          1.8020e+03, 3.9450e+03, 5.7190e+03, 2.2000e+02, 0.0000e+00],
#[Out]#         [1.0000e+00, 1.5000e+01, 7.5000e+01, 9.0000e+01, 1.5500e+02,
#[Out]#          2.8100e+02, 3.7700e+02, 5.2000e+01, 0.0000e+00, 0.0000e+00],
#[Out]#         [0.0000e+00, 3.0000e+00, 1.3000e+01, 1.7000e+01, 3.2000e+01,
#[Out]#          6.2000e+01, 1.0000e+01, 0.0000e+00, 0.0000e+00, 0.0000e+00],
#[Out]#         [0.0000e+00, 1.0000e+00, 0.0000e+00, 4.0000e+00, 1.9000e+01,
#[Out]#          7.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00]]),
#[Out]#  array([-2.227671 , -1.5341234, -0.8405758, -0.1470282,  0.5465194,
#[Out]#          1.240067 ,  1.9336146,  2.6271622,  3.3207098,  4.0142574,
#[Out]#          4.707805 ]),
#[Out]#  array([ 9.5907755 , 10.41923475, 11.247694  , 12.07615325, 12.9046125 ,
#[Out]#         13.73307175, 14.561531  , 15.38999025, 16.2184495 , 17.04690875,
#[Out]#         17.875368  ]),
#[Out]#  <matplotlib.collections.QuadMesh at 0x152abd6d4710>)
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'], sanders['phot_ks_mean_mag'], bins=[100,100])
#[Out]# (array([[0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         ...,
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.]]),
#[Out]#  array([-2.227671  , -2.15831624, -2.08896148, -2.01960672, -1.95025196,
#[Out]#         -1.8808972 , -1.81154244, -1.74218768, -1.67283292, -1.60347816,
#[Out]#         -1.5341234 , -1.46476864, -1.39541388, -1.32605912, -1.25670436,
#[Out]#         -1.1873496 , -1.11799484, -1.04864008, -0.97928532, -0.90993056,
#[Out]#         -0.8405758 , -0.77122104, -0.70186628, -0.63251152, -0.56315676,
#[Out]#         -0.493802  , -0.42444724, -0.35509248, -0.28573772, -0.21638296,
#[Out]#         -0.1470282 , -0.07767344, -0.00831868,  0.06103608,  0.13039084,
#[Out]#          0.1997456 ,  0.26910036,  0.33845512,  0.40780988,  0.47716464,
#[Out]#          0.5465194 ,  0.61587416,  0.68522892,  0.75458368,  0.82393844,
#[Out]#          0.8932932 ,  0.96264796,  1.03200272,  1.10135748,  1.17071224,
#[Out]#          1.240067  ,  1.30942176,  1.37877652,  1.44813128,  1.51748604,
#[Out]#          1.5868408 ,  1.65619556,  1.72555032,  1.79490508,  1.86425984,
#[Out]#          1.9336146 ,  2.00296936,  2.07232412,  2.14167888,  2.21103364,
#[Out]#          2.2803884 ,  2.34974316,  2.41909792,  2.48845268,  2.55780744,
#[Out]#          2.6271622 ,  2.69651696,  2.76587172,  2.83522648,  2.90458124,
#[Out]#          2.973936  ,  3.04329076,  3.11264552,  3.18200028,  3.25135504,
#[Out]#          3.3207098 ,  3.39006456,  3.45941932,  3.52877408,  3.59812884,
#[Out]#          3.6674836 ,  3.73683836,  3.80619312,  3.87554788,  3.94490264,
#[Out]#          4.0142574 ,  4.08361216,  4.15296692,  4.22232168,  4.29167644,
#[Out]#          4.3610312 ,  4.43038596,  4.49974072,  4.56909548,  4.63845024,
#[Out]#          4.707805  ]),
#[Out]#  array([ 9.5907755 ,  9.67362143,  9.75646735,  9.83931328,  9.9221592 ,
#[Out]#         10.00500513, 10.08785105, 10.17069698, 10.2535429 , 10.33638883,
#[Out]#         10.41923475, 10.50208068, 10.5849266 , 10.66777253, 10.75061845,
#[Out]#         10.83346438, 10.9163103 , 10.99915623, 11.08200215, 11.16484808,
#[Out]#         11.247694  , 11.33053993, 11.41338585, 11.49623178, 11.5790777 ,
#[Out]#         11.66192363, 11.74476955, 11.82761547, 11.9104614 , 11.99330732,
#[Out]#         12.07615325, 12.15899917, 12.2418451 , 12.32469102, 12.40753695,
#[Out]#         12.49038287, 12.5732288 , 12.65607472, 12.73892065, 12.82176658,
#[Out]#         12.9046125 , 12.98745842, 13.07030435, 13.15315028, 13.2359962 ,
#[Out]#         13.31884212, 13.40168805, 13.48453398, 13.5673799 , 13.65022582,
#[Out]#         13.73307175, 13.81591768, 13.8987636 , 13.98160952, 14.06445545,
#[Out]#         14.14730138, 14.2301473 , 14.31299323, 14.39583915, 14.47868508,
#[Out]#         14.561531  , 14.64437693, 14.72722285, 14.81006878, 14.8929147 ,
#[Out]#         14.97576063, 15.05860655, 15.14145248, 15.2242984 , 15.30714433,
#[Out]#         15.38999025, 15.47283618, 15.5556821 , 15.63852803, 15.72137395,
#[Out]#         15.80421988, 15.8870658 , 15.96991173, 16.05275765, 16.13560358,
#[Out]#         16.2184495 , 16.30129543, 16.38414135, 16.46698728, 16.5498332 ,
#[Out]#         16.63267913, 16.71552505, 16.79837098, 16.8812169 , 16.96406283,
#[Out]#         17.04690875, 17.12975468, 17.2126006 , 17.29544653, 17.37829245,
#[Out]#         17.46113838, 17.5439843 , 17.62683023, 17.70967615, 17.79252208,
#[Out]#         17.875368  ]),
#[Out]#  <matplotlib.collections.QuadMesh at 0x152abd67b650>)
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'], sanders['phot_ks_mean_mag'], bins=[100,100]);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=200);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[selection],
          sanders['phot_ks_mean_mag'][selection], bins=[100,100], vmax=500);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=[100,100], vmax=500);
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100], vmax=500);
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100], vmax=300);
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100], vmax=200);
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100], vmax=300);
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pmra'][highpm], sanders['pmdec'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -2) & significant & (sanders['pmra'] > -6) & (sanders['pmdec'] > -7)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pmra'][highpm], sanders['pmdec'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -4) & significant & (sanders['pmra'] > -6) & (sanders['pmdec'] > -9)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pmra'][highpm], sanders['pmdec'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -3) & significant & (sanders['pmra'] > -6) & (sanders['pmdec'] > -9)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pmra'][highpm], sanders['pmdec'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pmra']) > sanders['pmra_error'] * 2) &
                (np.abs(sanders['pmdec']) > sanders['pmdec_error'] * 2) &
                (sanders['pmra_error'] < 10) &
                (sanders['pmdec_error'] < 10)) |
               ((sanders['pmra_error'] < 1) &
                (sanders['pmdec_error'] < 1)))
highpm = (sanders['pmra'] < -2) & (sanders['pmdec'] < -3) & significant & (sanders['pmra'] > -6) & (sanders['pmdec'] > -9)
pl.scatter(sanders['pmra'][significant], sanders['pmdec'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pmra'][highpm], sanders['pmdec'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (RA) [mas yr$^{-1}$]")
pl.ylabel("PM (Dec) [mas yr$^{-1}$]");
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=[100,100], vmax=500);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100], vmax=300);
cm = pl.cm.viridis.copy()
cm.set_under(np.nan)
cm = pl.cm.viridis.copy()
cm.set_under((0,0,0,0))
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500, vmin=1);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500, vmin=1, cmap=cm);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500, vmin=1,
          norm=pl.cm.LogNorm(), cmap=cm);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500, vmin=1,
          norm=pl.colormaps.LogNorm(), cmap=cm);
from astropy.visualization import simple_norm
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100], vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=500, stretch='log'), cmap=cm);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=[100,100],# vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=500, stretch='log'), cmap=cm);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=[100,100],
          norm=simple_norm(data=[], vmin=1, vmax=500, stretch='log'), cmap=cm);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=[100,100],
          norm=simple_norm(data=[], vmin=1, vmax=500, stretch='log'), cmap=cm);
bins = [np.linspace(0.1,4), np.linspace(10, 17.5)]
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=bins,# vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=500, stretch='log'), cmap=cm);
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=bins,# vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=bins,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
from astropy.visualization import simple_norm
bins = [np.linspace(0.1,4,100), np.linspace(10, 17.5,100)]
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=bins,# vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1])
#[Out]# (17.5, 10.0)
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=bins,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1]);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=bins,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1]);
significant = (((np.abs(sanders['pml']) > sanders['pml_error'] * 2) &
                (np.abs(sanders['pmb']) > sanders['pmb_error'] * 2) &
                (sanders['pml_error'] < 10) &
                (sanders['pmb_error'] < 10)) |
               ((sanders['pml_error'] < 1) &
                (sanders['pmb_error'] < 1)))
highpm = (sanders['pml'] < -2) & (sanders['pmb'] < -3) & significant & (sanders['pml'] > -6) & (sanders['pmb'] > -9)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='olnge')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] < -2) & (sanders['pmb'] < -3) & significant & (sanders['pml'] > -6) & (sanders['pmb'] > -9)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='olnge')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] < -2) & (sanders['pmb'] < -3) & significant & (sanders['pml'] > -6) & (sanders['pmb'] > -9)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] < -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] > -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
cm = pl.cm.viridis.copy()
cm.set_under((0,0,0,0))
from astropy.visualization import simple_norm
bins = [np.linspace(0.1,4,100), np.linspace(10, 17.5,100)]
pl.hist2d(sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'],
          sanders['phot_ks_mean_mag'], bins=bins,# vmax=500, vmin=1,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1]);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[significant],
          sanders['phot_ks_mean_mag'][significant], bins=bins,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1]);
pl.hist2d((sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'])[highpm],
          sanders['phot_ks_mean_mag'][highpm], bins=bins,
          norm=simple_norm(data=[], vmin=1, vmax=1000, stretch='log'), cmap=cm);
pl.ylim(pl.ylim()[::-1]);
red_bright = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.0) & (sanders['phot_ks_mean_mag'] < 12)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][red_bright], sanders['pmb'][red_bright], s=1, alpha=0.1, c='orange')
pl.axis([-15, 15, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
red_bright = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.0) & (sanders['phot_ks_mean_mag'] < 12)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][red_bright], sanders['pmb'][red_bright], s=1, alpha=0.1, c='orange')
pl.axis([-15, 1, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
red_bright = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.0) & (sanders['phot_ks_mean_mag'] < 12)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.05)
pl.scatter(sanders['pml'][red_bright], sanders['pmb'][red_bright], s=1, alpha=0.1, c='orange')
pl.axis([-15, 1, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
red_bright = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.0) & (sanders['phot_ks_mean_mag'] < 12)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.05)
pl.scatter(sanders['pml'][red_bright], sanders['pmb'][red_bright], s=1, alpha=0.2, c='orange')
pl.axis([-15, 1, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
red_bright = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.0) & (sanders['phot_ks_mean_mag'] < 12)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.05)
pl.scatter(sanders['pml'][red_bright], sanders['pmb'][red_bright], s=1, alpha=0.5, c='orange')
pl.axis([-15, 1, -15, 15]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
np.linspace(1, -1)
#[Out]# array([ 1.        ,  0.95918367,  0.91836735,  0.87755102,  0.83673469,
#[Out]#         0.79591837,  0.75510204,  0.71428571,  0.67346939,  0.63265306,
#[Out]#         0.59183673,  0.55102041,  0.51020408,  0.46938776,  0.42857143,
#[Out]#         0.3877551 ,  0.34693878,  0.30612245,  0.26530612,  0.2244898 ,
#[Out]#         0.18367347,  0.14285714,  0.10204082,  0.06122449,  0.02040816,
#[Out]#        -0.02040816, -0.06122449, -0.10204082, -0.14285714, -0.18367347,
#[Out]#        -0.2244898 , -0.26530612, -0.30612245, -0.34693878, -0.3877551 ,
#[Out]#        -0.42857143, -0.46938776, -0.51020408, -0.55102041, -0.59183673,
#[Out]#        -0.63265306, -0.67346939, -0.71428571, -0.75510204, -0.79591837,
#[Out]#        -0.83673469, -0.87755102, -0.91836735, -0.95918367, -1.        ])
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(359.1, 358.42), np.linspace(-1/6, 1/6)])
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.42, 359.1, ), np.linspace(-1/6, 1/6)])
#[Out]# (array([[0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         ...,
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.],
#[Out]#         [0., 0., 0., ..., 0., 0., 0.]]),
#[Out]#  array([358.42      , 358.43387755, 358.4477551 , 358.46163265,
#[Out]#         358.4755102 , 358.48938776, 358.50326531, 358.51714286,
#[Out]#         358.53102041, 358.54489796, 358.55877551, 358.57265306,
#[Out]#         358.58653061, 358.60040816, 358.61428571, 358.62816327,
#[Out]#         358.64204082, 358.65591837, 358.66979592, 358.68367347,
#[Out]#         358.69755102, 358.71142857, 358.72530612, 358.73918367,
#[Out]#         358.75306122, 358.76693878, 358.78081633, 358.79469388,
#[Out]#         358.80857143, 358.82244898, 358.83632653, 358.85020408,
#[Out]#         358.86408163, 358.87795918, 358.89183673, 358.90571429,
#[Out]#         358.91959184, 358.93346939, 358.94734694, 358.96122449,
#[Out]#         358.97510204, 358.98897959, 359.00285714, 359.01673469,
#[Out]#         359.03061224, 359.0444898 , 359.05836735, 359.0722449 ,
#[Out]#         359.08612245, 359.1       ]),
#[Out]#  array([-0.16666667, -0.15986395, -0.15306122, -0.1462585 , -0.13945578,
#[Out]#         -0.13265306, -0.12585034, -0.11904762, -0.1122449 , -0.10544218,
#[Out]#         -0.09863946, -0.09183673, -0.08503401, -0.07823129, -0.07142857,
#[Out]#         -0.06462585, -0.05782313, -0.05102041, -0.04421769, -0.03741497,
#[Out]#         -0.03061224, -0.02380952, -0.0170068 , -0.01020408, -0.00340136,
#[Out]#          0.00340136,  0.01020408,  0.0170068 ,  0.02380952,  0.03061224,
#[Out]#          0.03741497,  0.04421769,  0.05102041,  0.05782313,  0.06462585,
#[Out]#          0.07142857,  0.07823129,  0.08503401,  0.09183673,  0.09863946,
#[Out]#          0.10544218,  0.1122449 ,  0.11904762,  0.12585034,  0.13265306,
#[Out]#          0.13945578,  0.1462585 ,  0.15306122,  0.15986395,  0.16666667]),
#[Out]#  <matplotlib.collections.QuadMesh at 0x152abb0ae480>)
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.42, 359.1, ), np.linspace(-1/6, 1/6)]);
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.6, 359.0, ), np.linspace(-0.1, 0.1)]);
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.6, 359.0, ), np.linspace(-0.1, 0.1)]);
pl.xlim(pl.xlim()[::-1]);
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.6, 359.0, ), np.linspace(-0.1, 0.1)]);
pl.xlim(pl.xlim()[::-1]);
pl.gca().set_aspect(1);
pl.hist2d(sanders['l'], sanders['b'], bins=[np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]);
pl.xlim(pl.xlim()[::-1]);
pl.gca().set_aspect(1);
pl.hist2d(sanders['l'][red_bright], sanders['b'][red_bright],
          bins=[np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]);
pl.xlim(pl.xlim()[::-1]);
pl.gca().set_aspect(1);
pl.hist2d(sanders['l'][highpm], sanders['b'][highpm],
          bins=[np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]);
pl.xlim(pl.xlim()[::-1]);
pl.gca().set_aspect(1);
glpth = '/orange/adamginsburg/cmz/glimpse_data/'
i1 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I1.fits')
i2 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I2.fits')
i2mi1 = i2 - i1
from astropy.io import fits
from astropy.wcs import WCS

glpth = '/orange/adamginsburg/cmz/glimpse_data/'
i1 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I1.fits')
i2 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I2.fits')
ww = wcs.WCS(fits.getheader(f'{glpth}/GLM_00000+0000_mosaic_I1.fits'))

i2mi1 = i2 - i1
from astropy.io import fits
from astropy.wcs import WCS

glpth = '/orange/adamginsburg/cmz/glimpse_data/'
i1 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I1.fits')
i2 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I2.fits')
ww = WCS(fits.getheader(f'{glpth}/GLM_00000+0000_mosaic_I1.fits'))

i2mi1 = i2 - i1
ww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : 0.0 0.0 
#[Out]# CRPIX : 4660.5 6760.5 
#[Out]# CD1_1 CD1_2  : -0.00033333333 0.0 
#[Out]# CD2_1 CD2_2  : 0.0 0.00033333333 
#[Out]# NAXIS : 9320  13520
import regions
from astropy import units as u
from astropy.coordinates import SkyCoord
coord = SkyCoord(358.7436209429, 0, unit=(u.deg, u.deg), frame='galactic')
l = 20*u.arcmin
w = 35*u.arcmin

coord_fova = SkyCoord('17:42:30.6943', '-30:10:46.49', unit=(u.hourangle, u.deg), frame='icrs')
coord_fovb = SkyCoord('17 42 48.5880', '-30 05 56.64', unit=(u.hourangle, u.deg), frame='icrs')
coord_fov_niriss = SkyCoord(coord_fovb.ra-4.*u.arcmin, coord_fovb.dec+5*u.arcmin, frame='icrs')
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)
msk = reg.to_pixel(ww).to_mask()
msk.cutout(i2mi1)
#[Out]# array([[  -4.6551   ,   -6.6052475,   -3.2984838, ...,   -6.7326517,
#[Out]#           -7.112303 ,   -5.259323 ],
#[Out]#        [ -15.529644 ,  -31.251339 ,  -14.847004 , ...,   -5.1043606,
#[Out]#           -4.4451323,   -6.6712914],
#[Out]#        [ -80.525734 , -124.15625  ,  -34.086395 , ...,   -2.859888 ,
#[Out]#           -3.2254915,   -5.633458 ],
#[Out]#        ...,
#[Out]#        [  -2.6722546,   -1.4483547,   -0.6702652, ...,   -2.9595985,
#[Out]#           -2.0076513,   -1.7916203],
#[Out]#        [  -6.4577255,   -2.469287 ,   -0.966589 , ...,   -2.7194052,
#[Out]#           -2.6180367,   -2.8759422],
#[Out]#        [  -4.777582 ,   -1.7507296,   -0.9599924, ...,   -3.361353 ,
#[Out]#           -3.4990234,   -1.9537239]], dtype=float32)
msk = reg.to_pixel(ww).to_mask()
msk.cutout(i2mi1)
msk.cutout(ww)
msk = reg.to_pixel(ww).to_mask()
msk.cutout(i2mi1)
msk.get_overlap_slices(ww)
msk = reg.to_pixel(ww).to_mask()
msk.cutout(i2mi1)
msk.get_overlap_slices(i1.shape)
#[Out]# ((slice(6259, 7261, None), slice(7554, 9305, None)),
#[Out]#  (slice(0, 1002, None), slice(0, 1751, None)))
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
pl.imshow( msk.cutout(i2mi1) )
#[Out]# <matplotlib.image.AxesImage at 0x152abab7b3b0>
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, min_percent=1))
#[Out]# <matplotlib.image.AxesImage at 0x152ab9494500>
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=95., min_percent=5))
#[Out]# <matplotlib.image.AxesImage at 0x152ab925ff80>
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99., min_percent=5))
pl.colorbar();
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99., min_percent=5))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='linear', max_percent=99.95, min_value=0))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='linear', max_percent=99.95, vmin=0))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='sqrt', max_percent=99.95, vmin=0))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=0))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
cutout[cutout<0] = 0
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=0))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
cutout[cutout<-5] = -5
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=-5))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
cutout[cutout<-25] = -25
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=-25))
pl.colorbar();
from astropy.io import fits
from astropy.wcs import WCS

glpth = '/orange/adamginsburg/cmz/glimpse_data/'
i1 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I1.fits')
i2 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I2.fits')
ww = WCS(fits.getheader(f'{glpth}/GLM_00000+0000_mosaic_I1.fits'))

i2mi1 = i2 - i1
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1)
cutout[cutout<-25] = -25
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=-25))
pl.colorbar();
from astropy.io import fits
from astropy.wcs import WCS

glpth = '/orange/adamginsburg/cmz/glimpse_data/'
i1 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I1.fits')
i2 = fits.getdata(f'{glpth}/GLM_00000+0000_mosaic_I2.fits')
ww = WCS(fits.getheader(f'{glpth}/GLM_00000+0000_mosaic_I1.fits'))

i2mi1 = i2 - i1
import regions
from astropy import units as u
from astropy.coordinates import SkyCoord
coord = SkyCoord(358.7436209429, 0, unit=(u.deg, u.deg), frame='galactic')
l = 20*u.arcmin
w = 35*u.arcmin

coord_fova = SkyCoord('17:42:30.6943', '-30:10:46.49', unit=(u.hourangle, u.deg), frame='icrs')
coord_fovb = SkyCoord('17 42 48.5880', '-30 05 56.64', unit=(u.hourangle, u.deg), frame='icrs')
coord_fov_niriss = SkyCoord(coord_fovb.ra-4.*u.arcmin, coord_fovb.dec+5*u.arcmin, frame='icrs')
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
cutout[cutout<-25] = -25
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=-25))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -10
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -8
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -5
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -1
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -2
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.figure(figsize=(12,8))
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -2
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.95, vmin=minval))
pl.colorbar();
msk = reg.to_pixel(ww).to_mask()
sl_big, sl_sml = msk.get_overlap_slices(i1.shape)
pl.figure(figsize=(12,8))
pl.subplot(projection=ww[sl_big])
cutout = msk.cutout(i2mi1).copy()
minval = -2
cutout[cutout<minval] = minval
pl.imshow(cutout, norm=simple_norm(cutout, stretch='log', max_percent=99.5, vmin=minval))
pl.colorbar();
for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=[np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
import scipy.ndimage as nd

cube_fn = '/orange/adamginsburg/galactic_plane_surveys/sedigism/sedigism/SgrE_13CO_smooth.fits'
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)

cube = SpectralCube.read(cube_fn).with_spectral_unit(u.km/u.s, velocity_convention='radio').subcube_from_regions([reg])
subcube = cube.spectral_slab(-220*u.km/u.s, -170*u.km/u.s)

noise = cube.spectral_slab(-300*u.km/u.s, -250*u.km/u.s).std(axis=0)
peak = subcube.max(axis=0)
spatial_mask = peak > 3*noise
subcube_masked = subcube.with_mask(spatial_mask)
mom0_masked = subcube_masked.moment0()
mom0 = subcube.moment0()
resampled = nd.gaussian_filter(mom0.value, 1)

def percentile_list(data, percents=[91]):
    return [np.nanpercentile(data, p) for p in percents]

plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='white')
import scipy.ndimage as nd
from spectral_cube import SpectralCube

cube_fn = '/orange/adamginsburg/galactic_plane_surveys/sedigism/sedigism/SgrE_13CO_smooth.fits'
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)

cube = SpectralCube.read(cube_fn).with_spectral_unit(u.km/u.s, velocity_convention='radio').subcube_from_regions([reg])
subcube = cube.spectral_slab(-220*u.km/u.s, -170*u.km/u.s)

noise = cube.spectral_slab(-300*u.km/u.s, -250*u.km/u.s).std(axis=0)
peak = subcube.max(axis=0)
spatial_mask = peak > 3*noise
subcube_masked = subcube.with_mask(spatial_mask)
mom0_masked = subcube_masked.moment0()
mom0 = subcube.moment0()
resampled = nd.gaussian_filter(mom0.value, 1)

def percentile_list(data, percents=[91]):
    return [np.nanpercentile(data, p) for p in percents]

plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='white')
import scipy.ndimage as nd
from spectral_cube import SpectralCube
import pylab as plt

cube_fn = '/orange/adamginsburg/galactic_plane_surveys/sedigism/sedigism/SgrE_13CO_smooth.fits'
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)

cube = SpectralCube.read(cube_fn).with_spectral_unit(u.km/u.s, velocity_convention='radio').subcube_from_regions([reg])
subcube = cube.spectral_slab(-220*u.km/u.s, -170*u.km/u.s)

noise = cube.spectral_slab(-300*u.km/u.s, -250*u.km/u.s).std(axis=0)
peak = subcube.max(axis=0)
spatial_mask = peak > 3*noise
subcube_masked = subcube.with_mask(spatial_mask)
mom0_masked = subcube_masked.moment0()
mom0 = subcube.moment0()
resampled = nd.gaussian_filter(mom0.value, 1)

def percentile_list(data, percents=[91]):
    return [np.nanpercentile(data, p) for p in percents]

plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='white')
#[Out]# <matplotlib.contour.QuadContourSet at 0x152ab3f94d70>
pl.imshow(resample)
pl.imshow(resampled)
#[Out]# <matplotlib.image.AxesImage at 0x152ab07fe5d0>
import scipy.ndimage as nd
from spectral_cube import SpectralCube
import pylab as plt

cube_fn = '/orange/adamginsburg/galactic_plane_surveys/sedigism/sedigism/SgrE_13CO_smooth.fits'
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)

cube = SpectralCube.read(cube_fn).with_spectral_unit(u.km/u.s, velocity_convention='radio').subcube_from_regions([reg])
subcube = cube.spectral_slab(-220*u.km/u.s, -170*u.km/u.s)

noise = cube.spectral_slab(-300*u.km/u.s, -250*u.km/u.s).std(axis=0)
peak = subcube.max(axis=0)
spatial_mask = peak > 3*noise
subcube_masked = subcube.with_mask(spatial_mask)
mom0_masked = subcube_masked.moment0()
mom0 = subcube.moment0()
resampled = nd.gaussian_filter(mom0.value, 1)

def percentile_list(data, percents=[91]):
    return [np.nanpercentile(data, p) for p in percents]

plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='lk')
import scipy.ndimage as nd
from spectral_cube import SpectralCube
import pylab as plt

cube_fn = '/orange/adamginsburg/galactic_plane_surveys/sedigism/sedigism/SgrE_13CO_smooth.fits'
reg = regions.RectangleSkyRegion(center=coord, width=w, height=l)

cube = SpectralCube.read(cube_fn).with_spectral_unit(u.km/u.s, velocity_convention='radio').subcube_from_regions([reg])
subcube = cube.spectral_slab(-220*u.km/u.s, -170*u.km/u.s)

noise = cube.spectral_slab(-300*u.km/u.s, -250*u.km/u.s).std(axis=0)
peak = subcube.max(axis=0)
spatial_mask = peak > 3*noise
subcube_masked = subcube.with_mask(spatial_mask)
mom0_masked = subcube_masked.moment0()
mom0 = subcube.moment0()
resampled = nd.gaussian_filter(mom0.value, 1)

def percentile_list(data, percents=[91]):
    return [np.nanpercentile(data, p) for p in percents]

plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='k')
#[Out]# <matplotlib.contour.QuadContourSet at 0x152ab0d75610>
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.crpix = [0,0]
histww.crval = [358.6, -0.1]
histww.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.crpix = [0,0]
histww.crval = [358.6, -0.1]
histww.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    #pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.crpix = [0,0]
histww.crval = [358.6, -0.1]
histww.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : '' '' 
#[Out]# CRVAL : 0.0 0.0 
#[Out]# CRPIX : 0.0 0.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 1.0 1.0 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0,0]
histww.wcs.crval = [358.6, -0.1]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : 358.6 -0.1 
#[Out]# CRPIX : 0.0 0.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.004040404040381418 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0,0]
histww.wcs.crval = [-1.4, -0.1]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 -0.1 
#[Out]# CRPIX : 0.0 0.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.004040404040381418 0.0020202020202020193 
#[Out]# NAXIS : 0  0
ww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : 0.0 0.0 
#[Out]# CRPIX : 4660.5 6760.5 
#[Out]# CD1_1 CD1_2  : -0.00033333333 0.0 
#[Out]# CD2_1 CD2_2  : 0.0 0.00033333333 
#[Out]# NAXIS : 9320  13520
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.004040404040381418 0.0020202020202020193 
#[Out]# NAXIS : 0  0
np.diff(np.linspace(358.6, 359.0, 100))
#[Out]# array([0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404, 0.0040404,
#[Out]#        0.0040404, 0.0040404, 0.0040404])
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    pl.hist2d(sanders['l'][sel], sanders['b'][sel],
              bins=bins, transform=histww);
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
np.histogram2d(sanders['l'][sel], sanders['b'][sel],
          bins=bins);
np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
#[Out]# (array([[ 1.,  5.,  2., ...,  6.,  6., 10.],
#[Out]#         [ 3.,  4.,  3., ...,  7.,  4.,  3.],
#[Out]#         [ 4.,  2.,  3., ...,  7.,  6.,  5.],
#[Out]#         ...,
#[Out]#         [ 0.,  0.,  1., ...,  0.,  0.,  0.],
#[Out]#         [ 0.,  0.,  1., ...,  0.,  0.,  0.],
#[Out]#         [ 0.,  1.,  2., ...,  0.,  0.,  0.]]),
#[Out]#  array([358.6       , 358.6040404 , 358.60808081, 358.61212121,
#[Out]#         358.61616162, 358.62020202, 358.62424242, 358.62828283,
#[Out]#         358.63232323, 358.63636364, 358.64040404, 358.64444444,
#[Out]#         358.64848485, 358.65252525, 358.65656566, 358.66060606,
#[Out]#         358.66464646, 358.66868687, 358.67272727, 358.67676768,
#[Out]#         358.68080808, 358.68484848, 358.68888889, 358.69292929,
#[Out]#         358.6969697 , 358.7010101 , 358.70505051, 358.70909091,
#[Out]#         358.71313131, 358.71717172, 358.72121212, 358.72525253,
#[Out]#         358.72929293, 358.73333333, 358.73737374, 358.74141414,
#[Out]#         358.74545455, 358.74949495, 358.75353535, 358.75757576,
#[Out]#         358.76161616, 358.76565657, 358.76969697, 358.77373737,
#[Out]#         358.77777778, 358.78181818, 358.78585859, 358.78989899,
#[Out]#         358.79393939, 358.7979798 , 358.8020202 , 358.80606061,
#[Out]#         358.81010101, 358.81414141, 358.81818182, 358.82222222,
#[Out]#         358.82626263, 358.83030303, 358.83434343, 358.83838384,
#[Out]#         358.84242424, 358.84646465, 358.85050505, 358.85454545,
#[Out]#         358.85858586, 358.86262626, 358.86666667, 358.87070707,
#[Out]#         358.87474747, 358.87878788, 358.88282828, 358.88686869,
#[Out]#         358.89090909, 358.89494949, 358.8989899 , 358.9030303 ,
#[Out]#         358.90707071, 358.91111111, 358.91515152, 358.91919192,
#[Out]#         358.92323232, 358.92727273, 358.93131313, 358.93535354,
#[Out]#         358.93939394, 358.94343434, 358.94747475, 358.95151515,
#[Out]#         358.95555556, 358.95959596, 358.96363636, 358.96767677,
#[Out]#         358.97171717, 358.97575758, 358.97979798, 358.98383838,
#[Out]#         358.98787879, 358.99191919, 358.9959596 , 359.        ]),
#[Out]#  array([-0.1       , -0.0979798 , -0.0959596 , -0.09393939, -0.09191919,
#[Out]#         -0.08989899, -0.08787879, -0.08585859, -0.08383838, -0.08181818,
#[Out]#         -0.07979798, -0.07777778, -0.07575758, -0.07373737, -0.07171717,
#[Out]#         -0.06969697, -0.06767677, -0.06565657, -0.06363636, -0.06161616,
#[Out]#         -0.05959596, -0.05757576, -0.05555556, -0.05353535, -0.05151515,
#[Out]#         -0.04949495, -0.04747475, -0.04545455, -0.04343434, -0.04141414,
#[Out]#         -0.03939394, -0.03737374, -0.03535354, -0.03333333, -0.03131313,
#[Out]#         -0.02929293, -0.02727273, -0.02525253, -0.02323232, -0.02121212,
#[Out]#         -0.01919192, -0.01717172, -0.01515152, -0.01313131, -0.01111111,
#[Out]#         -0.00909091, -0.00707071, -0.00505051, -0.0030303 , -0.0010101 ,
#[Out]#          0.0010101 ,  0.0030303 ,  0.00505051,  0.00707071,  0.00909091,
#[Out]#          0.01111111,  0.01313131,  0.01515152,  0.01717172,  0.01919192,
#[Out]#          0.02121212,  0.02323232,  0.02525253,  0.02727273,  0.02929293,
#[Out]#          0.03131313,  0.03333333,  0.03535354,  0.03737374,  0.03939394,
#[Out]#          0.04141414,  0.04343434,  0.04545455,  0.04747475,  0.04949495,
#[Out]#          0.05151515,  0.05353535,  0.05555556,  0.05757576,  0.05959596,
#[Out]#          0.06161616,  0.06363636,  0.06565657,  0.06767677,  0.06969697,
#[Out]#          0.07171717,  0.07373737,  0.07575758,  0.07777778,  0.07979798,
#[Out]#          0.08181818,  0.08383838,  0.08585859,  0.08787879,  0.08989899,
#[Out]#          0.09191919,  0.09393939,  0.0959596 ,  0.0979798 ,  0.1       ]))
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.004040404040381418 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
[bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
#[Out]# [0.004040404040381418, 0.0020202020202020193]
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.004040404040381418 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.0020100502512718776 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H.T)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.0020100502512718776 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H.T)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w', transform=mom0.wcs)
    break
histww
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    pl.imshow(H.T)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w', transform=ax.get_transform(mom0.wcs))
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.0020100502512718776 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w', transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
    break
histww
#[Out]# WCS Keywords
#[Out]# 
#[Out]# Number of WCS axes: 2
#[Out]# CTYPE : 'GLON-CAR' 'GLAT-CAR' 
#[Out]# CRVAL : -1.4 0.0 
#[Out]# CRPIX : 0.0 50.0 
#[Out]# PC1_1 PC1_2  : 1.0 0.0 
#[Out]# PC2_1 PC2_2  : 0.0 1.0 
#[Out]# CDELT : 0.0020100502512718776 0.0020202020202020193 
#[Out]# NAXIS : 0  0
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(1, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii} < H-K < {ii+0.1}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w', transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(0.5, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < H-K < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-15, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    sel =  (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < H-K < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-15, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    sel =  (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < H-K < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-15, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    sel =  (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < $\delta\ell$ < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-15, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    sel =  (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < $\\delta\\ell$ < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.01)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.01, c='orange')
pl.axis([-15, 15, -2.5, 2.5]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.1)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.1, c='orange')
pl.axis([-1, 15, -2.5, 2.5]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.1)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.1, c='orange')
pl.axis([-15, 1, -2.5, 2.5]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.1)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.1, c='orange')
pl.axis([-15, 1, -3.5, 3.5]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
significant = (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1)))
highpm = (sanders['pml'] > -10) & (sanders['pmb'] < 2) & significant & (sanders['pml'] < -4) & (sanders['pmb'] > -2)
pl.scatter(sanders['pml'][significant], sanders['pmb'][significant], s=1, alpha=0.1)
pl.scatter(sanders['pml'][highpm], sanders['pmb'][highpm], s=1, alpha=0.1, c='orange')
pl.axis([-15, 1, -4.5, 4.5]);
pl.xlabel("PM (l) [mas yr$^{-1}$]")
pl.ylabel("PM (b) [mas yr$^{-1}$]");
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 50)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-10, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.5) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < 1.8)
    sel &= (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < $\\delta\\ell$ < {ii+0.1:0.1f}, 1.5 < H-K < 1.8")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
bins = [np.linspace(358.6, 359.0, 100), np.linspace(-0.1, 0.1, 50)]

histww = WCS()
histww.wcs.crpix = [0, len(bins[1])/2]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-10, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.5) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < 1.8)
    sel &= (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < $\\delta\\ell$ < {ii+0.1:0.1f}, 1.5 < H-K < 1.8")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
get_ipython().run_line_magic('mkdir', 'figures')
bins = [np.linspace(358.6, 359.0, 50), np.linspace(-0.1, 0.1, 25)]

histww = WCS()
histww.wcs.crpix = [0, len(bins[1])/2]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(-10, -1, 1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > 1.5) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < 1.8)
    sel &= (((np.abs(sanders['pml']) > sanders['epml'] * 2) &
                (np.abs(sanders['pmb']) > sanders['epmb'] * 2) &
                (sanders['epml'] < 10) &
                (sanders['epmb'] < 10)) |
               ((sanders['epml'] < 1) &
                (sanders['epmb'] < 1))) & (sanders['pml'] > ii) & (sanders['pml'] < ii + 1) 
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < $\\delta\\ell$ < {ii+1:0.1f}, 1.5 < H-K < 1.8")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
    pl.savefig(f'figures/histogram_dl_{ii:d}to{ii+1:d}_1.5ltHKlt1.8.png', bbox_inches='tight')
bins = [np.linspace(358.6, 359.0, 200), np.linspace(-0.1, 0.1, 100)]

histww = WCS()
histww.wcs.crpix = [0, 50]
histww.wcs.crval = [-1.4, 0]
histww.wcs.cdelt = [bins[0][1] - bins[0][0], bins[1][1] - bins[1][0]]
histww.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']

for ii in np.arange(0.5, 3, 0.1):
    pl.figure()
    sel = (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] > ii) & (sanders['phot_h_mean_mag']-sanders['phot_ks_mean_mag'] < ii + 0.1)
    ax = pl.subplot(projection=histww)
    H, x, y = np.histogram2d(sanders['l'][sel], sanders['b'][sel], bins=bins)
    ax.imshow(H.T)
    pl.title(f"{ii:0.1f} < H-K < {ii+0.1:0.1f}")
    pl.xlim(pl.xlim()[::-1]);
    #pl.gca().set_aspect(1);
    axis = pl.axis()
    plt.contour(resampled, levels=[percentile_list(resampled, [80])[0]], colors='w',
                linewidths=[0.5],
                transform=ax.get_transform(mom0.wcs))
    pl.axis(axis)
    pl.savefig(f'figures/histogram_{ii:0.1f}ltHKlt{ii+0.1:0.1f}.png', bbox_inches='tight')