notebook.community

Edit and run



In [7]:

    
import matplotlib.pyplot as plt
import numpy as np
from copy import deepcopy
import pickle
from astropy.coordinates import SkyCoord
from dustmaps.config import config
config['data_dir'] = '/home/weissj3/Desktop/MWTools/Data'

import dustmaps.sfd
dustmaps.sfd.fetch()
import astropy.units as u









    



Downloading SFD data file to /home/weissj3/Desktop/MWTools/Data/sfd/SFD_dust_4096_ngp.fits
Checking existing file to see if MD5 sum matches ...
File exists. Not overwriting.
Downloading SFD data file to /home/weissj3/Desktop/MWTools/Data/sfd/SFD_dust_4096_sgp.fits
Checking existing file to see if MD5 sum matches ...
File exists. Not overwriting.



In [9]:

    
filenames = ["N12toN10", "N15toN12", "N17toN15", "N20toN17", "P10toP11", "P11toP12", "P12toP13", "P13top14", "P14top15", "P15toP16", "P16toP17", "P17toP20"]
starLen = 0
for name in filenames:
    data = {"ID":[], "g":[], "r":[], "g_err":[], "r_err":[], "l":[], "b":[]}
    f = open("../../SummerREUDiskWork/public_html/PanStarrsMSTO/PS_MSTO/" + name + ".csv", "r")
    f.readline()


    for line in f:
        ln = list(map(float, line.split(",")))
        data["ID"].append(ln[0])
        data["g"].append(ln[1])
        data["r"].append(ln[2])
        data["g_err"].append(ln[3])
        data["r_err"].append(ln[4])
        data["l"].append(ln[5])
        data["b"].append(ln[6])
    print("Stars In File %s: %d" % (name, len(data['ID'])))
    coords = SkyCoord(data['l'],data['b'],  unit='deg', frame='galactic')

    sfd = dustmaps.sfd.SFDQuery()
    extinctionValues = sfd(coords) #3.303 is conversion from E(B-v) to sdss g
    data['g_0'] = np.array(data['g']) - (extinctionValues * 3.172)
    data['r_0'] = np.array(data['r']) - (extinctionValues * 2.271)
    data['(g-r)_0'] = np.array(data['g_0']) - np.array(data['r_0'])
    
    f.close()


    mask = np.where((np.array(data["(g-r)_0"]) > 0.1) & (np.array(data["(g-r)_0"]) < 0.3)) 

    data["ID"] = np.array(data["ID"])[mask]
    data["g"] = np.array(data["g"])[mask]
    data["r"] = np.array(data["r"])[mask]
    data["g_err"] = np.array(data["g_err"])[mask]
    data["r_err"] = np.array(data["r_err"])[mask]
    data["l"] = np.array(data["l"])[mask]
    data["b"] = np.array(data["b"])[mask]
    data['g_0'] = np.array(data['g_0'])[mask]
    data['r_0'] = np.array(data['r_0'])[mask]
    data['(g-r)_0'] = np.array(data['(g-r)_0'])[mask]
    
    a = open("../../SummerREUDiskWork/public_html/PanStarrsMSTO/ColorCut/" + name + ".pk", "wb")
    pickle.dump(data, a)
    a.close()
    del data









    



Stars In File N12toN10: 10259083
Stars In File N15toN12: 2126541
Stars In File N17toN15: -5937026
Stars In File N20toN17: 1919483
Stars In File P10toP11: -3338635
Stars In File P11toP12: -375503
Stars In File P12toP13: -415586
Stars In File P13top14: -237238
Stars In File P14top15: -276338
Stars In File P15toP16: -215792
Stars In File P16toP17: -213184






    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-9-496988ef2501> in <module>()
     18     print("Stars In File %s: %d" % (name, len(data['ID']) - starLen))
     19     starLen = len(data['ID'])
---> 20     coords = SkyCoord(data['l'],data['b'],  unit='deg', frame='galactic')
     21 
     22     sfd = dustmaps.sfd.SFDQuery()

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/sky_coordinate.py in __init__(self, copy, *args, **kwargs)
    223         # creating the internal self._sky_coord_frame object
    224         args = list(args)  # Make it mutable
--> 225         kwargs = self._parse_inputs(args, kwargs)
    226 
    227         frame = kwargs['frame']

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/sky_coordinate.py in _parse_inputs(self, args, kwargs)
    434                                                                       repr_attr_names, units):
    435                     attr_class = frame.representation.attr_classes[repr_attr_name]
--> 436                     coord_kwargs[frame_attr_name] = attr_class(arg, unit=unit)
    437 
    438             else:

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/angles.py in __new__(cls, angle, unit, **kwargs)
    504         if isinstance(angle, Longitude):
    505             raise TypeError("A Latitude angle cannot be created from a Longitude angle")
--> 506         self = super().__new__(cls, angle, unit=unit, **kwargs)
    507         self._validate_angles()
    508         return self

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/angles.py in __new__(cls, angle, unit, dtype, copy)
    105                   not (isinstance(angle, np.ndarray) and
    106                        angle.dtype.kind not in 'SUVO')):
--> 107                 angle = [Angle(x, unit, copy=False) for x in angle]
    108 
    109         return super().__new__(cls, angle, unit, dtype=dtype, copy=copy)

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/angles.py in <listcomp>(.0)
    105                   not (isinstance(angle, np.ndarray) and
    106                        angle.dtype.kind not in 'SUVO')):
--> 107                 angle = [Angle(x, unit, copy=False) for x in angle]
    108 
    109         return super().__new__(cls, angle, unit, dtype=dtype, copy=copy)

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/angles.py in __new__(cls, angle, unit, dtype, copy)
    107                 angle = [Angle(x, unit, copy=False) for x in angle]
    108 
--> 109         return super().__new__(cls, angle, unit, dtype=dtype, copy=copy)
    110 
    111     @staticmethod

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/quantity.py in __new__(cls, value, unit, dtype, copy, order, subok, ndmin)
    384 
    385         value = value.view(cls)
--> 386         value._set_unit(value_unit)
    387         if unit is value_unit:
    388             return value

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/coordinates/angles.py in _set_unit(self, unit)
    129 
    130     def _set_unit(self, unit):
--> 131         super()._set_unit(self._convert_unit_to_angle_unit(unit))
    132 
    133     @property

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/quantity.py in _set_unit(self, unit)
   1715 
   1716     def _set_unit(self, unit):
-> 1717         if unit is None or not unit.is_equivalent(self._equivalent_unit):
   1718             raise UnitTypeError(
   1719                 "{0} instances require units equivalent to '{1}'"

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/core.py in is_equivalent(self, other, equivalencies)
    787         other = Unit(other, parse_strict='silent')
    788 
--> 789         return self._is_equivalent(other, equivalencies)
    790 
    791     def _is_equivalent(self, other, equivalencies=[]):

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/core.py in _is_equivalent(self, other, equivalencies)
    799 
    800         if (self._get_physical_type_id() ==
--> 801                 other._get_physical_type_id()):
    802             return True
    803         elif len(equivalencies):

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/core.py in _get_physical_type_id(self)
    535         """
    536         unit = self.decompose()
--> 537         r = zip([x.name for x in unit.bases], unit.powers)
    538         # bases and powers are already sorted in a unique way
    539         # r.sort()

/usr/local/lib/python3.4/dist-packages/astropy-3.1.dev21527-py3.4-linux-x86_64.egg/astropy/units/core.py in powers(self)
    587         Return the powers of the unit.
    588         """
--> 589         return [1]
    590 
    591     def to_string(self, format=unit_format.Generic):

KeyboardInterrupt:



In [ ]:

    
#Read in Pannstars Data from file (JEFF'S DATA (NEW))
#Don't use this if you can avoid it. Eats all memory
data = {"ID":[], "g":[], "r":[], "g_err":[], "r_err":[], "l":[], "b":[]}
filenames = ["N12toN10.csv", "N15toN12.csv", "N17toN15.csv", "N20toN17.csv", "P10toP11.csv", "P11toP12.csv", "P12toP13.csv", "P13top14.csv", "P14top15.csv", "P15toP16.csv", "P16toP17.csv", "P17toP20.csv"]
starLen = 0
for name in filenames:
    f = open("../../SummerREUDiskWork/public_html/PanStarrsMSTO/PS_MSTO/" + name, "r")
    f.readline()


    for line in f:
        ln = list(map(float, line.split(",")))
        data["ID"].append(ln[0])
        data["g"].append(ln[1])
        data["r"].append(ln[2])
        data["g_err"].append(ln[3])
        data["r_err"].append(ln[4])
        data["l"].append(ln[5])
        data["b"].append(ln[6])
    print("Stars In File %s: %d" % (name, len(data['ID']) - starLen))
    starLen = len(data['ID'])

    f.close()

print("Read in %d stars" % len(data["ID"]))









    



Stars In File N12toN10.csv: 10259083
Stars In File N15toN12.csv: 12385624
Stars In File N17toN15.csv: 6448598
Stars In File N20toN17.csv: 8368081
Stars In File P10toP11.csv: 5029446
Stars In File P11toP12.csv: 4653943
Stars In File P12toP13.csv: 4238357
Stars In File P13top14.csv: 4001119
Stars In File P14top15.csv: 3724781
Stars In File P15toP16.csv: 3508989



In [ ]:

    
#Read in Pannstars Data from file (SARAH'S DATA (OLD))
data = {"ID":[], "g":[], "r":[], "g_err":[], "r_err":[], "l":[], "b":[]}
filenames = ["25L100-20B-10.csv", "25L10010B20.csv", "100L210-20B-10.csv", "100L21010B20.csv", "210L330-20B-10.csv", "210L33010B20.csv", "330L25allB.csv"]
starLen = 0
for name in filenames:
    f = open("../../SummerREUDiskWork/public_html/SarahGonzalez2018/PANSTARRS DR1 Data/Data with Extinctions and Applied Color Cut/" + name, "r")
    f.readline()


    for line in f:
        ln = list(map(float, line.split(",")))
        data["ID"].append(ln[1])
        data["g"].append(ln[2])
        data["r"].append(ln[3])
        data["g_err"].append(ln[4])
        data["r_err"].append(ln[5])
        data["l"].append(ln[6])
        data["b"].append(ln[7])
    print("Stars In File %s: %d" % (name, len(data['ID']) - starLen))
    starLen = len(data['ID'])

    f.close()

print("Read in %d stars" % len(data["ID"]))



In [3]:

    
import gc

gc.collect()









    Out[3]:





0



In [5]:

    
for i in np.arange(len(data["l"])):
    data["l"][i] = float(data["l"][i] > 180) * -(360 - data["l"][i]) +  float(data["l"][i] < 180) * (data["l"][i])



In [6]:

    
#Plot Pannstars data to make sure it was read in correctly
#Looks like there is a lot of dust
plt.figure(1, figsize=(16, 12))
mask = np.where((np.array(data["g"]) > 16.0) & (np.array(data["g"]) < 22.5)) 

plt.plot(np.array(data["l"])[mask],np.array(data["b"])[mask], 'o', ms=0.1, alpha=0.02)
plt.show()



In [6]:

    
#Slice the disk into l bins

for i in np.arange(0, 360, 2.5):
    temp = []
    mask = np.where((np.array(data["l"]) > i) & (np.array(data["l"]) < i+2.5)) 
    temp.append(np.array(data["l"])[mask])
    temp.append(np.array(data["b"])[mask])
    temp.append(np.array(data["g"])[mask])
    a = open("../Data/SlicedDisk/l_cuts/l-%d.pickle" % int((i+1.25)*100), "wb")
    pickle.dump(temp, a)
    a.close()



In [50]:

    
#Plot to check cuts
plt.plot(np.array(l_cut)[3, 0],np.array(l_cut)[3, 1], 'o', ms=0.1, alpha=0.5)
plt.show()



In [20]:

    
#Cut into b bins
for i in np.arange(0, 360, 2.5):
    a = open("../Data/SlicedDisk/l_cuts/l-%d.pickle" % int((i+1.25)*100), "rb")
    loadedData = pickle.load(a)
    a.close()
    for j in np.arange(-20, 20, 2.5):
        temp = []
        mask = np.where((np.array(loadedData[1]) > j) & (np.array(loadedData[1]) < j+2.5)) 
        temp.append(np.array(loadedData[0])[mask])
        temp.append(np.array(loadedData[1])[mask])
        temp.append(np.array(loadedData[2])[mask])
        b = open("../Data/SlicedDisk/lb_cuts/l-%d-b-%d.pickle" % (int((i+1.25)*100), int((j+1.25)*100)), "wb")
        pickle.dump(temp, b)
        b.close()
    a.close()



In [25]:

    
#Plot distance histogram
b = open("../Data/SlicedDisk/lb_cuts/l-1125-b-1625.pickle", "rb")
temp = pickle.load(b)
b.close()


plt.hist(temp[2])
plt.show()



In [35]:

    
#Convert pickle files to histogram files for optimizer
for i in np.arange(0, 360, 2.5):
    for j in np.arange(-20, 20, 2.5):
        a = open("../Data/SlicedDisk/lb_cuts/l-%d-b-%d.pickle" % (int((i+1.25)*100), int((j+1.25)*100)), "rb")
        loadedData = pickle.load(a)
        a.close()
        hist = np.histogram(loadedData[2], bins=13, range=(16.0, 22.5))
        b = open("../Data/SlicedDisk/lb_histograms/l-%d-b-%d.hist" % (int((i+1.25)*100), int((j+1.25)*100)), "w")
        b.write('0')
        for k in hist[0]:
            b.write(', ' + str(k))
        b.close()
    a.close()



In [ ]:

    
#Plot Results
for i in np.arange(0, 360, 2.5):
    for j in np.arange(-20, 20, 2.5):
        a = open("../Data/SlicedDisk/lb_cuts/l-%d-b-%d.pickle" % (int((i+1.25)*100), int((j+1.25)*100)), "rb")
        loadedData = pickle.load(a)
        a.close()
        hist = np.histogram(loadedData[2], bins=13, range=(16.0, 22.5))
        b = open("../Data/SlicedDisk/lb_histograms/l-%d-b-%d.hist" % (int((i+1.25)*100), int((j+1.25)*100)), "w")
        b.write('0')
        for k in hist[0]:
            b.write(', ' + str(k))
        b.close()
    a.close()