How to Make your own OM10 sim input catalogs in 4 easy steps


In [1]:
%matplotlib inline
from __future__ import division
import om10,os
import numpy as np
import matplotlib.pyplot as plt
import triangle

1. Load the OM10 pre-simulated lenses

Load the catalog of simulated lenses. Pandas interface between FITS tables and DataFrames only exists in the dev version of astropy, so use native astropy tables for now!


In [2]:
db = om10.DB(catalog='$OM10_DIR/data/qso_mock.fits')
print db.lenses


  LENSID  FLAGTYPE NIMG ZLENS VELDISP  ... SIGCRIT  DSLUM  L_I REFF REFF_T
--------- -------- ---- ----- -------- ... ------- ------- --- ---- ------
    14428        0    2 0.062 264.2409 ...     0.0 35956.6 0.0  0.0    0.0
    21703        0    2 0.072 165.1962 ...     0.0  8193.9 0.0  0.0    0.0
    31711        0    2 0.082 155.9553 ...     0.0 35956.6 0.0  0.0    0.0
    34331        0    2 0.084 172.5838 ...     0.0 14498.4 0.0  0.0    0.0
    51810        0    2 0.096 244.9063 ...     0.0 26239.2 0.0  0.0    0.0
    83414        0    2 0.114   154.17 ...     0.0 11904.1 0.0  0.0    0.0
   115074        0    2 0.128 113.5011 ...     0.0 7101.16 0.0  0.0    0.0
   125027        0    2  0.13 247.1724 ...     0.0 16302.5 0.0  0.0    0.0
   147709        0    2 0.138 200.9093 ...     0.0 17748.4 0.0  0.0    0.0
   180815        0    2 0.148 185.3532 ...     0.0 14029.4 0.0  0.0    0.0
   202550        0    4 0.154  181.134 ...     0.0 8794.67 0.0  0.0    0.0
      ...      ...  ...   ...      ... ...     ...     ... ...  ...    ...
113167635        0    2 1.936 250.6109 ...     0.0 47637.1 0.0  0.0    0.0
115614728        0    4  1.96 242.1029 ...     0.0 38304.5 0.0  0.0    0.0
118666840        0    2  1.99 207.9697 ...     0.0 27265.9 0.0  0.0    0.0
120546345        0    2 2.008 281.8383 ...     0.0 28090.1 0.0  0.0    0.0
123598414        0    2 2.038 199.5262 ...     0.0 31000.1 0.0  0.0    0.0
129002676        0    2  2.09 236.0478 ...     0.0 39270.7 0.0  0.0    0.0
138366386        0    2  2.18 227.5097 ...     0.0 36169.7 0.0  0.0    0.0
157075130        0    2 2.358 258.8213 ...     0.0 29955.8 0.0  0.0    0.0
168229053        0    2 2.464 200.9093 ...     0.0 48954.0 0.0  0.0    0.0
185790887        0    2  2.63 215.7744 ...     0.0 42184.9 0.0  0.0    0.0
193208710        0    2   2.7 259.4179 ...     0.0 37235.8 0.0  0.0    0.0

2. Place lenses on sky

Read in positions, redshifts, and magnitudes of CFHTLS galaxies from file, then match with simulated lenses by nearest neighbor.


In [16]:
db.get_sky_positions(input_cat='$OM10_DIR/data/CFHTLS_LRGs.txt')
db.assign_sky_positions()


om10.DB: read in LRG sky position data from  /Users/mbaumer/pybin/OM10/data/CFHTLS_LRGs.txt
Mean LRG RA,DEC,z,i =  34.3807307408 -7.09643979181 0.612498 21.63331883
om10.DB: number of LRGs stored =  10000
15658

In [17]:
print db.lenses


  LENSID  FLAGTYPE NIMG ZLENS VELDISP  ... L_I REFF REFF_T   RA   DEC 
--------- -------- ---- ----- -------- ... --- ---- ------ ----- -----
    14428        0    2 0.062 264.2409 ... 0.0  0.0    0.0 -99.0 -99.0
    21703        0    2 0.072 165.1962 ... 0.0  0.0    0.0 -99.0 -99.0
    31711        0    2 0.082 155.9553 ... 0.0  0.0    0.0 -99.0 -99.0
    34331        0    2 0.084 172.5838 ... 0.0  0.0    0.0 -99.0 -99.0
    51810        0    2 0.096 244.9063 ... 0.0  0.0    0.0 -99.0 -99.0
    83414        0    2 0.114   154.17 ... 0.0  0.0    0.0 -99.0 -99.0
   115074        0    2 0.128 113.5011 ... 0.0  0.0    0.0 -99.0 -99.0
   125027        0    2  0.13 247.1724 ... 0.0  0.0    0.0 -99.0 -99.0
   147709        0    2 0.138 200.9093 ... 0.0  0.0    0.0 -99.0 -99.0
   180815        0    2 0.148 185.3532 ... 0.0  0.0    0.0 -99.0 -99.0
   202550        0    4 0.154  181.134 ... 0.0  0.0    0.0 -99.0 -99.0
      ...      ...  ...   ...      ... ... ...  ...    ...   ...   ...
113167635        0    2 1.936 250.6109 ... 0.0  0.0    0.0 -99.0 -99.0
115614728        0    4  1.96 242.1029 ... 0.0  0.0    0.0 -99.0 -99.0
118666840        0    2  1.99 207.9697 ... 0.0  0.0    0.0 -99.0 -99.0
120546345        0    2 2.008 281.8383 ... 0.0  0.0    0.0 -99.0 -99.0
123598414        0    2 2.038 199.5262 ... 0.0  0.0    0.0 -99.0 -99.0
129002676        0    2  2.09 236.0478 ... 0.0  0.0    0.0 -99.0 -99.0
138366386        0    2  2.18 227.5097 ... 0.0  0.0    0.0 -99.0 -99.0
157075130        0    2 2.358 258.8213 ... 0.0  0.0    0.0 -99.0 -99.0
168229053        0    2 2.464 200.9093 ... 0.0  0.0    0.0 -99.0 -99.0
185790887        0    2  2.63 215.7744 ... 0.0  0.0    0.0 -99.0 -99.0
193208710        0    2   2.7 259.4179 ... 0.0  0.0    0.0 -99.0 -99.0

~1/3 of the objects have an RA,DEC of (-99,-99). This is because there were no galaxies found in the reference sample with matching magnitudes/redshifts...


In [18]:
print len(db.lenses[db.lenses['RA']> 0]), len(db.lenses)


10482 15658

3. Simulate photometry

Add photometry by finding similar objects in SDSS/WISE. This step takes a while, so one can also load up a pre-painted catalog in this directory.


In [3]:
db.paint(lrg_input_cat='$OM10_DIR/data/LRGo.txt',qso_input_cat='$OM10_DIR/data/QSOo.txt')


[[  1.49700100e-03   3.80608700e+02   6.37357900e+00 ...,   1.54550000e+01
    1.26080000e+01   8.39200000e+00]
 [  1.67441600e-03   1.21834900e+02   4.75030300e-03 ...,   1.43040000e+01
    1.20260000e+01   8.69300000e+00]
 [  2.50279000e-03   1.15952900e+02   5.66436100e-01 ...,   1.50150000e+01
    1.18490000e+01   8.09800000e+00]
 ..., 
 [  9.68880800e-01   2.76625200e+02   3.75018900e-01 ...,   1.22730000e+01
    9.37900000e+00   7.61400000e+00]
 [  9.70397100e-01   1.88763100e+02   2.22591200e-01 ...,   1.04410000e+01
    9.34700000e+00   8.04300000e+00]
 [  9.74084200e-01   1.34568700e+02   1.12254600e+00 ...,   1.44720000e+01
    1.14880000e+01   7.81100000e+00]]

In [5]:
print db.lenses['MAGG_SRC','MAGR_SRC','MAGI_SRC','MAGZ_SRC', \
        'MAGW1_SRC','MAGW2_SRC','MAGW3_SRC','MAGW4_SRC', 'SDSS_FLAG_SRC']


MAGG_SRC MAGR_SRC MAGI_SRC MAGZ_SRC ... MAGW3_SRC MAGW4_SRC SDSS_FLAG_SRC
-------- -------- -------- -------- ... --------- --------- -------------
21.52536 19.80611 19.39259 19.39907 ...    11.469     8.302           0.0
18.79825 18.53064 18.58717 18.68878 ...    10.772     7.631           0.0
21.52536 19.80611 19.39259 19.39907 ...    11.469     8.302           0.0
18.82549 18.89963  18.7185 18.68868 ...    11.547     7.988           0.0
19.57941 19.07229 18.80365 18.75532 ...    10.684     7.736           0.0
19.83415 19.75743 19.55235 19.52748 ...    11.984     8.757           0.0
17.77648 17.56068 17.56235 17.57965 ...     9.534     7.051           0.0
18.89332 19.02151 18.81176 18.63794 ...    12.025     7.831           0.0
18.73156 18.67111 18.64709 18.36551 ...    10.662     7.907           0.0
19.58352 19.60018 19.39213 19.38316 ...    12.274     9.069           0.0
18.91262 18.72567 18.79178 18.80928 ...    11.093     8.727           0.0
     ...      ...      ...      ... ...       ...       ...           ...
22.24258 21.39198 20.47871 20.05188 ...    10.934     8.377           0.0
20.58433 18.81212 18.12769  17.7807 ...    10.409     8.275           0.0
 19.0098 18.76121 18.72146 18.77421 ...    11.302     8.278           0.0
19.88507  19.4411 19.25687 19.35262 ...    11.824     8.515           0.0
 18.5316 17.84844   17.828 17.73992 ...    10.633      8.27           0.0
21.38347 20.61231 20.29492 19.78409 ...    11.798     8.417           0.0
20.71608 19.15028  18.9268 18.92575 ...    11.527     8.349           0.0
19.39724 18.53825  18.4894 18.44785 ...    11.188     8.057           0.0
20.70699 19.58138 18.97668 19.06188 ...    10.472     7.084           0.0
19.77423 19.57165 19.11099  18.9079 ...     10.55     8.705           0.0
21.46115  19.8115 19.41338 19.45662 ...    11.904     8.884           0.0

In [21]:
from scipy.stats import itemfreq
itemfreq(db.lenses['SDSS_FLAG_SRC'])


Out[21]:
array([[  0.00000000e+00,   1.55570000e+04],
       [  1.00000000e+00,   1.01000000e+02]])

With 2-param matching, we had 1941 matched and 13717 unmatched. Matching in z alone, we have 15557 matched and 101 unmatched!


In [6]:
#write out cat
db.lenses.write('my_redshift_matching.fits',format='fits')

In [20]:
#pre-load painted, positioned catalog for this NB
db = om10.DB(catalog='$OM10_DIR/notebooks/painted_positioned_catalog.fits')

4. Select the sample you want


In [9]:
db.select_random(maglim=21.4,area=30000.0,IQ=1.0)


om10.DB: selection yields  995  lenses

5. Output flat catalog


In [10]:
out_cat = db.make_sim_input_catalog()

In [27]:
print out_cat.colnames


['LENSID', 'RA', 'DEC', 'XIMG', 'YIMG', 'G', 'R', 'I', 'Z']

In [29]:
out_cat


Out[29]:
LENSIDRADEC...IZ
24079482.034.05474379-6.450411733...12.508565076412.183990741
24079482.034.0549478456-6.45047631633...4.214924370814.11585783407
24079482.034.0546904844-6.45040751078...5.489539142255.39047260551
6656185.034.67705221-4.991044321...19.344063776919.0341610445
6656185.034.6773356822-4.991022571...18.417436088918.3693331904
6656185.034.6769906267-4.99101537656...19.213065498519.1649626
20752770.034.39178125-4.567873844...11.798738561911.3421201723
20752770.034.3915580556-4.56807381622...16.89908926316.7743128604
20752770.034.3914792778-4.56792806622...16.816495510616.691719108
20752770.034.3920063889-4.56787626067...18.895919110918.7711427083
20752770.034.3915629444-4.567640344...17.660613985517.5358375829
..................
15074352.034.0659280011-5.45992719244...3.682938078663.63268666809
15074352.034.0653177233-5.45980335911...4.740128909364.68987749879
1706588.034.42226001-9.280916479...18.466151596118.5783651985
1706588.034.4224300378-9.28081942344...3.439113211543.37951839149
1706588.034.4221307878-9.28098206233...2.083214749972.02361992992
75322.0-99.0-99.0...16.936760061419.8230577549
75322.0-99.0002504167-99.0000246944...19.070111501818.8130474862
75322.0-98.9999406944-99.0000049722...20.724845149620.467781134
14303170.034.02372569-10.05982628...19.082850956518.0463669018
14303170.034.0238728289-10.0599386967...18.161006357218.1765026382
14303170.034.0233704956-10.05952178...17.973101075217.9885973562

In [14]:
positioned_cat = out_cat[out_cat['RA']>0] #let's only take galaxies that were matched to a location

In [27]:
positioned_cat.write('sim_input_sample_cat.fits',format='fits')

Exploring our positioned, painted catalog


In [13]:
for var in ['G','R','I','Z']:
    plt.figure()
    plt.hist(positioned_cat[var][positioned_cat['XIMG'] == 0],bins=20,histtype='step',color='Orange',label='Lenses')
    plt.hist(positioned_cat[var][positioned_cat['XIMG'] != 0],bins=20,histtype='step',color='Cyan',label='Sources')
    plt.legend()
    plt.xlabel(var)



In [25]:
plt.scatter(positioned_cat['RA'],positioned_cat['DEC'],color='Orange')
plt.title('Assigned Locations of Simulated Lenses in CFHT footprint (?)')
plt.xlabel('RA')
plt.ylabel('DEC')


Out[25]:
<matplotlib.text.Text at 0x10e09c050>

In [40]:
plt.hist(db.lenses['APMAG_I'])


Out[40]:
(array([   17.,    58.,   225.,   652.,  1632.,  3109.,  4203.,  3746.,
         1739.,   277.]),
 array([ 13.8852615 ,  15.01993596,  16.15461042,  17.28928488,
         18.42395934,  19.5586338 ,  20.69330826,  21.82798272,
         22.96265718,  24.09733164,  25.2320061 ]),
 <a list of 10 Patch objects>)

In [47]:
plt.scatter(db.lenses['ZLENS'],db.lenses['VELDISP'])
plt.xlabel('z')
plt.ylabel('Sigma')


Out[47]:
<matplotlib.text.Text at 0x11142d410>

In [10]:
from astropy.table import Table
sloan = Table.read('../data/LRGo.txt',format='ascii')
print sloan


   col1      col2       col3       col4    ... col11  col12  col13  col14
---------- -------- ----------- ---------- ... ------ ------ ------ -----
0.08777109 101.4772    29.66101  0.5684491 ...  14.19 13.869 10.882 7.935
0.08915185 119.6037    9.038554   7.443291 ... 13.661 13.513 10.027 8.289
0.08907133 147.3923    7.934169   6.422607 ...   13.4 13.309 10.887 8.071
0.07815439 126.4102    7.423387   9.857728 ... 14.067 13.945 10.701  8.28
0.07330353 88.20443   0.0426188   11.35288 ... 13.257 13.123  9.269 7.541
0.06905825 111.7736    0.478357   0.925322 ... 13.279 13.146 10.043 7.802
0.07845496 101.5035    29.65839    2.14771 ... 14.225  14.16 10.509 8.675
0.09718453 97.15182 0.001023908    6.28618 ... 13.667 13.506 12.078 8.031
0.08647919 140.8284    2.346406   6.159338 ... 13.301 13.532  9.953 7.732
0.09598833 83.64037    1.026456   8.836071 ...  13.39 13.367 10.215 7.662
0.04160196 92.87505 0.001052819   2.211381 ... 13.141 13.133 11.277 7.552
       ...      ...         ...        ... ...    ...    ...    ...   ...
 0.9410933 232.3272    1.177355   1.185546 ... 14.993 14.984 11.613 8.275
  0.908578 361.5846    3.397827   1.297545 ... 14.432 14.289 11.946 8.292
 0.9398239  186.413    1.173412   1.034775 ... 14.484 14.249 12.301  8.17
 0.9688808 276.6252   0.3750189 0.06943174 ...  13.67 12.273  9.379 7.614
 0.9250459 317.8982     2.56627   3.663085 ... 15.438 15.274 11.357 7.778
 0.9403334   390.15    5.761228   0.722985 ... 14.477 14.237 11.264 7.886
   0.96223 221.9585    0.882306   1.054062 ... 14.746 14.703 11.807 8.461
    0.9353  383.818   0.9183709   2.107691 ... 15.056 15.369 12.052 8.558
 0.9703971 188.7631   0.2225912 0.03742412 ... 11.406 10.441  9.347 8.043
 0.9446794  177.801   0.4517192   1.468057 ... 13.982 13.854 12.673 8.791
 0.9740842 134.5687    1.122546   2.336176 ... 14.548 14.472 11.488 7.811

In [11]:
plt.scatter(db.sample['ZLENS'],db.sample['VELDISP'],color='Orange')
plt.scatter(sloan['col1'],sloan['col2'],color='Blue')
plt.xlabel('z')
plt.ylabel('Sigma')


Out[11]:
<matplotlib.text.Text at 0x10aa81550>

In [ ]: