Merge pull request #45 from LSSTDESC/u/jchiang/instcat_bug_fixes

bug-fixes and added flag to disable proper motion

scripts/generateDc2InstCat.py → scripts/protoDC2/generateDc2InstCat.py

@@ -3,6 +3,7 @@
 import os
 import numpy as np
 import gzip
+import h5py
 
 from lsst.sims.catUtils.exampleCatalogDefinitions import PhoSimCatalogPoint
 from lsst.sims.catalogs.definitions import InstanceCatalog
@@ -12,20 +13,41 @@
 
 class MaskedPhoSimCatalogPoint(PhoSimCatalogPoint):
 
+    disable_proper_motion = False
+
     min_mag = None
 
     column_outputs = ['prefix', 'uniqueId', 'raPhoSim', 'decPhoSim', 'maskedMagNorm', 'sedFilepath',
-                      'redshift', 'gamma1', 'gamma2', 'kappa', 'raOffset', 'decOffset',
+                      'redshift', 'shear1', 'shear2', 'kappa', 'raOffset', 'decOffset',
                       'spatialmodel', 'internalExtinctionModel',
                       'galacticExtinctionModel', 'galacticAv', 'galacticRv']
 
+    protoDc2_half_size = 2.5*np.pi/180.
+
     @cached
     def get_maskedMagNorm(self):
         raw_norm = self.column_by_name('phoSimMagNorm')
         if self.min_mag is None:
             return raw_norm
         return np.where(raw_norm<self.min_mag, self.min_mag, raw_norm)
 
+    @cached
+    def get_inProtoDc2(self):
+        ra_values = self.column_by_name('raPhoSim')
+        ra = np.where(ra_values < np.pi, ra_values, ra_values - 2.*np.pi)
+        dec = self.column_by_name('decPhoSim')
+        return np.where((ra > -self.protoDc2_half_size) &
+                        (ra < self.protoDc2_half_size) &
+                        (dec > -self.protoDc2_half_size) &
+                        (dec < self.protoDc2_half_size), 1, None)
+
+    def column_by_name(self, colname):
+        if (self.disable_proper_motion and
+            (colname.startswith('properMotion')
+             or colname == 'radialVelocity')):
+            return np.zeros(len(self.column_by_name('raJ2000')), dtype=np.float)
+        return super(MaskedPhoSimCatalogPoint, self).column_by_name(colname)
+
 
 class BrightStarCatalog(PhoSimCatalogPoint):
 
@@ -52,21 +74,21 @@ def get_isBright(self):
     parser.add_argument('--id', type=int, nargs='+',
                         default=None,
                         help='obsHistID to generate InstanceCatalog for (a list)')
-    parser.add_argument('--dither', type=str,
-                        default='True',
-                        help='whether or not to apply dithering (true/false; default true)')
+    parser.add_argument('--disable_dithering', default=False,
+                        action='store_true',
+                        help='flag to disable dithering')
     parser.add_argument('--min_mag', type=float, default=10.0,
                         help='the minimum magintude for stars')
     parser.add_argument('--fov', type=float, default=2.0,
                         help='field of view radius in degrees')
+    parser.add_argument('--disable_proper_motion', default=False,
+                        action='store_true',
+                        help='flag to disable proper motion')
     args = parser.parse_args()
 
     obshistid_list = args.id
     opsimdb = args.db
     out_dir = args.out
-    dither_switch = True
-    if args.dither.lower()[0] == 'f':
-        dither_switch = False
 
     from lsst.sims.catUtils.utils import ObservationMetaDataGenerator
 
@@ -98,8 +120,7 @@ def get_isBright(self):
                                                         boundLength=args.fov)
 
         obs = obs_list[0]
-        if dither_switch:
-            print 'dithering'
+        if not args.disable_dithering:
             obs.pointingRA = np.degrees(obs.OpsimMetaData['randomDitherFieldPerVisitRA'])
             obs.pointingDec = np.degrees(obs.OpsimMetaData['randomDitherFieldPerVisitDec'])
             rotSky = _getRotSkyPos(obs._pointingRA, obs._pointingDec, obs,
@@ -121,10 +142,12 @@ def get_isBright(self):
             output.write('includeobj %s.gz\n' % gal_name)
             #output.write('includeobj %s.gz\n' % agn_name)
 
-        star_cat = MaskedPhoSimCatalogPoint(star_db, obs_metadata=obs)
+        star_cat = MaskedPhoSimCatalogPoint(star_db, obs_metadata=obs,
+                                            cannot_be_null=['inProtoDc2'])
         star_cat.phoSimHeaderMap = phosim_header_map
         bright_cat = BrightStarCatalog(star_db, obs_metadata=obs, cannot_be_null=['isBright'])
         star_cat.min_mag = args.min_mag
+        star_cat.disable_proper_motion = args.disable_proper_motion
         bright_cat.min_mag = args.min_mag
 
         from lsst.sims.catalogs.definitions import parallelCatalogWriter
@@ -143,9 +166,9 @@ def get_isBright(self):
         cat.write_catalog(os.path.join(out_dir, gal_name), chunk_size=100000,
                           write_mode='a', write_header=False)
 
-        for orig_name in (star_name, gal_name, agn_name):
+        for orig_name in (star_name, gal_name):
             full_name = os.path.join(out_dir, orig_name)
             with open(full_name, 'r') as input_file:
-                with gzip.open(full_name+'.gz', 'w') as output_file:
+                with gzip.open(full_name+'.gz', 'wb') as output_file:
                     output_file.writelines(input_file)
             os.unlink(full_name)

scripts/protoDC2/partition_visits.py

@@ -0,0 +1,11 @@
+import pandas as pd
+
+df = pd.read_pickle('protoDC2_visits.pkl')
+
+for filt in 'ugrizy':
+    df_filt = df.query('filter=="%s"' % filt)
+    df_filt.sort_values('expMJD')
+    with open('protoDC2_visits_%s-band.txt' % filt, 'w') as output:
+        for irow in df_filt.index:
+            output.write('%i  %s\n' % (df_filt.loc[irow]['obsHistID'],
+                                       df_filt.loc[irow]['expMJD']))

scripts/protoDC2/plot_visits.py

@@ -0,0 +1,41 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import gzip
+
+plt.ion()
+
+def plot_galaxies(catalog, limit=1000000, alpha=0.01, color='grey'):
+    ra, dec = [], []
+    with gzip.open(catalog) as input_:
+        for i, line in enumerate(input_):
+            if i > limit:
+                break
+            tokens = line.split()
+            ra_val = float(tokens[2])
+            if ra_val > 180:
+                ra_val -= 360.
+            ra.append(ra_val)
+            dec.append(float(tokens[3]))
+    plt.errorbar(ra, dec, fmt='.', alpha=alpha, color=color)
+
+def plot_visits(df, xname, yname, side=2.5, ms=2, alpha=1):
+    for filt in 'ugrizy':
+        band = df.query('filter=="%s"' % filt)
+        label = '%s band, %i visits' % (filt, len(band))
+        plt.errorbar(band[xname], band[yname], fmt='.', label=label, ms=ms,
+                     alpha=alpha)
+    plt.xlabel(xname)
+    plt.ylabel(yname)
+    plt.legend(fontsize='x-small', loc=0)
+    plt.plot((-side, -side, side, side, -side),
+             (-side, side, side, -side, -side))
+
+df = pd.read_pickle('protoDC2_visits.pkl')
+
+plt.figure()
+plot_galaxies('gal_cat_138143.txt.gz')
+plot_visits(df, 'randomDitherFieldPerVisitRA', 'randomDitherFieldPerVisitDec')
+
+plt.title('minion_1016_sqlite_new_dithers.db, %i total visits' % len(df))
+plt.savefig('protoDC2_visits.png')

scripts/protoDC2/protoDC2_valid_obsids.py

@@ -0,0 +1,66 @@
+import pickle
+import numpy as np
+import pandas as pd
+from lsst.sims.catUtils.utils import ObservationMetaDataGenerator
+
+def fov_overlaps_protoDC2(x, y, half_size=2.5, radius=1.77):
+    if x > half_size and y > half_size:
+        return (x - half_size)**2 + (y - half_size)**2 < radius**2
+    elif x < -half_size and y > half_size:
+        return (x + half_size)**2 + (y - half_size)**2 < radius**2
+    elif x < -half_size and y < -half_size:
+        return (x + half_size)**2 + (y + half_size)**2 < radius**2
+    elif x > half_size and y < -half_size:
+        return (x - half_size)**2 + (y + half_size)**2 < radius**2
+    outer_box_size = half_size + radius
+    return (-outer_box_size < x and x < outer_box_size and
+            -outer_box_size < y and y < outer_box_size)
+
+opsim_db = 'minion_1016_sqlite_new_dithers.db'
+obs_gen = ObservationMetaDataGenerator(opsim_db)
+
+# Inclusive fov is 2.1 deg radius; protoDC2 is 5x5 box centered on
+# (RA, Dec) = (0, 0), with vertical sides aligned N-S and horizontal
+# sides aligned E-W.  Use 20x20 box to ensure all viable dithered
+# visits are considered
+dxy = 10
+fov = 2.1
+radius = fov
+obs_list = obs_gen.getObservationMetaData(fieldRA=(0, dxy),
+                                          fieldDec=(-dxy, dxy),
+                                          boundLength=radius)
+obs_list.extend(obs_gen.getObservationMetaData(fieldRA=(360-dxy, 360.),
+                                               fieldDec=(-dxy, dxy),
+                                               boundLength=radius))
+
+df = pd.DataFrame(columns=['obsHistID', 'fieldRA', 'fieldDec',
+                           'randomDitherFieldPerVisitRA',
+                           'randomDitherFieldPerVisitDec',
+                           'filter', 'fieldID', 'propID', 'expMJD'])
+
+for obs in obs_list:
+    ditheredRA = (obs.summary['OpsimMetaData']['randomDitherFieldPerVisitRA']
+                  *180./np.pi)
+    if ditheredRA > 180:
+        ditheredRA -= 360.
+    ditheredDec = (obs.summary['OpsimMetaData']['randomDitherFieldPerVisitDec']
+                   *180./np.pi)
+
+    if not fov_overlaps_protoDC2(ditheredRA, ditheredDec):
+        continue
+
+    fieldRA = obs.summary['OpsimMetaData']['fieldRA']*180./np.pi
+    if fieldRA > 180:
+        fieldRA -= 360.
+    df.loc[len(df)] = (obs.summary['OpsimMetaData']['obsHistID'],
+                       fieldRA,
+                       obs.summary['OpsimMetaData']['fieldDec']*180./np.pi,
+                       ditheredRA, ditheredDec,
+                       obs.summary['OpsimMetaData']['filter'],
+                       obs.summary['OpsimMetaData']['fieldID'],
+                       obs.summary['OpsimMetaData']['propID'],
+                       obs.summary['OpsimMetaData']['expMJD'])
+
+pickle.dump(df, open('protoDC2_visits.pkl', 'wb'), protocol=2)
+
+print(len(df))