A demonstration of sims_movingObjects

This notebook is also available in the sims_movingObjects github repository, in the examples directory: https://github.com/lsst/sims_movingObjects/blob/master/examples/Demo%20sims_movingObjects.ipynb

[1]:

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

[2]:

import lsst.sims.movingObjects as mo
import lsst.sims.maf.db as db
from lsst.sims.maf.batches import getColMap

For most typical use of sims_movingObjects, you just need to use the script ‘makeLSSTobs.py’ :

(lsst-scipipe) [ewok:~] lynnej% makeLSSTobs.py --help

usage: makeLSSTobs.py [-h] [--opsimDb OPSIMDB] [--orbitFile ORBITFILE]
                      [--outDir OUTDIR] [--obsFile OBSFILE]
                      [--sqlConstraint SQLCONSTRAINT]
                      [--obsMetadata OBSMETADATA] [--footprint FOOTPRINT]
                      [--rFov RFOV] [--xTol XTOL] [--yTol YTOL]
                      [--roughTol ROUGHTOL] [--obsType OBSTYPE]
                      [--obsCode OBSCODE] [--tStep TSTEP] [--ephMode EPHMODE]
                      [--prelimEphMode PRELIMEPHMODE] [--ephType EPHTYPE]

Generate moving object detections.

optional arguments:
  -h, --help            show this help message and exit
  --opsimDb OPSIMDB     Opsim output db file (example: kraken_2026.db).
                        Default None.
  --orbitFile ORBITFILE
                        File containing the moving object orbits. See
                        https://github.com/lsst/oorb/blob/lsst-
                        dev/python/README.rst for additional documentation on
                        the orbit file format. Default None.
  --outDir OUTDIR       Output directory for moving object detections. Default
                        '.'
  --obsFile OBSFILE     Output file name for moving object observations.
                        Default will build outDir/opsimRun_orbitFile_obs.txt.
  --sqlConstraint SQLCONSTRAINT
                        SQL constraint to use to select data from opsimDb.
                        Default no constraint.
  --obsMetadata OBSMETADATA
                        Additional metadata to write into output file. The
                        default metadata will combine the opsimDb name, the
                        sqlconstraint, and the name of the orbit file;
                        obsMetadata is an optional addition.
  --footprint FOOTPRINT
                        Type of footprint to use to identify observations of
                        each object. Options are 'circle', 'rectangle', or
                        'camera' (apply camera footprint). Default is 'circle'
                        (which will then have a 1.75 deg radius).
  --rFov RFOV           If using a circular footprint, this is the radius of
                        the FOV (in degrees). Default 1.75 degrees.
  --xTol XTOL           If using a rectangular footprint, this is the
                        tolerance in the RA direction (in degrees). Default is
                        5 degrees.
  --yTol YTOL           If using a rectangular footprint, this is the
                        tolerance in the Dec direction (in degrees). Default
                        is 3 degrees.
  --roughTol ROUGHTOL   If using direct/exact ephemeris generation, this is
                        the tolerance for the preliminary matches between
                        ephemerides and pointings (in degrees). Default 20
                        degrees.
  --obsType OBSTYPE     Method for generating observations: 'direct' or
                        'linear'. Linear will use linear interpolation between
                        a grid of ephemeris points. Direct will first generate
                        rough ephemerides, look for observations within
                        roughTol of these points, and then generate exact
                        ephemerides at those times. Default 'direct'.
  --obsCode OBSCODE     Observatory code for generating observations. Default
                        is I11 (Cerro Pachon).
  --tStep TSTEP         Timestep between ephemeris generation for either the
                        first (rough) stage of direct ephemeris generation or
                        the grid for linear interpolation ephemerides. Default
                        1 day.
  --ephMode EPHMODE     2body or nbody mode for ephemeris generation. Default
                        is nbody.
  --prelimEphMode PRELIMEPHMODE
                        Use either 2body or nbody for preliminary ephemeris
                        generation in the rough stage for DirectObs. Default
                        2body.
  --ephType EPHTYPE     Generate either 'basic' or 'full' ephemerides from
                        OOrb. See https://github.com/lsst/oorb/blob/lsst-
                        dev/python/README.rst for detailsof the contents of
                        'full' or 'basic' ephemerides. Default basic.

This runs the steps below, as a whole, to produce an output file with the observations.

But let’s break down a little bit of what’s happening behind the scenes.

Reading the orbit file.

[3]:

# First the orbit file - this is read into an Orbits object.
orbitDir = '.'
orbitfile = 'neos_100.s3m'

orbits = mo.Orbits()
orbits.readOrbits(os.path.join(orbitDir, orbitfile))

[4]:

# The actual orbits are checked to see if they contain required columns:
for orbit_type in orbits.dataCols:
    print(orbit_type, orbits.dataCols[orbit_type])

COM ['objId', 'q', 'e', 'inc', 'Omega', 'argPeri', 'tPeri', 'epoch', 'H', 'g', 'sed_filename']
KEP ['objId', 'a', 'e', 'inc', 'Omega', 'argPeri', 'meanAnomaly', 'epoch', 'H', 'g', 'sed_filename']
CART ['objId', 'x', 'y', 'z', 'xdot', 'ydot', 'zdot', 'epoch', 'H', 'g', 'sed_filename']

[5]:

# and if 'standard but not quite the same' orbit column names are used,
# the column names are swapped to these standard names.

# Note that columns 'objId' = anything (including strings), and columns H / g / sed_filename are not required.
# If H is not specified, the value of 20 will be used
# If g is not specified, the value of 0.15 will be used
# If sed_filename is not specified (it should match a filename in $SIMS_MOVINGOBJECTS_DIR/data), one will be assigned.

help(orbits.assignSed)

Help on method assignSed in module lsst.sims.movingObjects.orbits:

assignSed(orbits, randomSeed=None) method of lsst.sims.movingObjects.orbits.Orbits instance
    Assign either a C or S type SED, depending on the semi-major axis of the object.
    P(C type) = 0 (a<2); 0.5*a - 1 (2<a<4); 1 (a > 4),
    based on figure 23 from Ivezic et al 2001 (AJ, 122, 2749).

    Parameters
    ----------
    orbits : pandas.DataFrame, pandas.Series or numpy.ndarray
       Array-like object containing orbital parameter information.

    Returns
    -------
    numpy.ndarray
        Array containing the SED type for each object in 'orbits'.


[6]:

# The final orbit data is kept in a class attribute: orbits:
# note that additional data in the orbit file (such as otype) is retained.
orbits.orbits[0:5]

[6]:
objId otype q e inc Omega argPeri tPeri H g epoch model sed_filename
0 S1000000 S3M 3.01822 0.05208 22.56035 211.00286 335.42134 57801.64493 14.20 0.15 57231.0 V150728 C.dat
1 S1000001 S3M 2.10974 0.07518 4.91571 209.40298 322.66447 56722.78779 20.57 0.15 57231.0 V150728 S.dat
2 S1000002 S3M 2.80523 0.07777 1.24945 112.52284 139.86858 56406.44304 14.65 0.15 57231.0 V150728 S.dat
3 S1000003 S3M 2.10917 0.13219 1.46615 266.54621 232.24412 56980.11676 19.58 0.15 57231.0 V150728 S.dat
4 S1000004 S3M 2.17676 0.19949 12.92422 162.14580 192.22312 56808.90252 10.56 0.15 57231.0 V150728 C.dat

Connecting to (and reading from) the opsim data file.

[7]:

# Look at the MAF documentation for more information on the lsst.sims.maf.db.OpsimDatabase class and methods.
opsimDir = '/Users/lynnej/opsim/db'
opsimrun = 'baseline2018b'
opsdb = db.OpsimDatabase(os.path.join(opsimDir, opsimrun + '.db'))
colmap = getColMap(opsdb)

[8]:

# Basically, we will just read all of these requested columns from the database: (for now we'll stick to night < 100)
reqcols = [colmap['mjd'], colmap['night'], colmap['ra'], colmap['dec'],
               'rotSkyPos', colmap['filter'], colmap['exptime'], colmap['seeingEff'],
               colmap['seeingGeom'], colmap['fiveSigmaDepth'], 'solarElong']
degreesIn = colmap['raDecDeg']
print(reqcols)

['observationStartMJD', 'night', 'fieldRA', 'fieldDec', 'rotSkyPos', 'filter', 'visitExposureTime', 'seeingFwhmEff', 'seeingFwhmGeom', 'fiveSigmaDepth', 'solarElong']

[9]:

simdata = opsdb.fetchMetricData(reqcols, sqlconstraint='night < 100')

[10]:

# Simdata is a numpy recarray - let's just use pandas to make it easier to read here.
pd.DataFrame(simdata[0:3])

[10]:
observationStartMJD night fieldRA fieldDec rotSkyPos filter visitExposureTime seeingFwhmEff seeingFwhmGeom fiveSigmaDepth solarElong
0 59853.016794 1 305.088793 -24.889283 180.000000 z 30.0 0.692029 0.620848 23.348066 113.652908
1 59853.018021 1 304.770694 35.966846 180.786087 z 30.0 1.182104 1.023690 22.308947 113.957998
2 59853.020185 1 305.482200 -62.802603 181.134222 z 30.0 0.766303 0.681901 23.178827 99.684211

Generating observations.

There is more than one way to actually generate observations, but let’s look at the ‘direct’ approach. This is implemented in the class “DirectObs”, which provides higher level methods that basically do the following:

  • generate (rough) ephemerides between the start and end of the survey data, on stepsizes of tStep (typically, 1 day)
  • identify pointings which come within roughTol of these rough ephemerides
  • generate exact ephemerides at the times of those pointings, and test if they fall within ‘footprint’
  • ‘footprint’ can be circle, rectangle or camera – ‘circle’ just checks if an object falls within rFov of the center of the pointing, ‘rectangle’ checks if RA x cos(dec) is within xTol and Dec within yTol of the pointing center, and camera applies the actual LSST camera footprint to check if an object falls on active silicon.
  • for the pointings where the object was observed, then calculate trailing losses and write the output to disk (including a color term, appropriate for the filter of each observation).

(Note: there’s plenty of room to improve on the procedure above, which could both speed things up and improve accuracy .. if this is something you’re interested in contributing to, please consider writing a new Obs class, preferably inheriting from the BaseObs class).

[11]:

# Using the "DirectObs" class.
rFov = 1.75  # in degrees
obscode = 'I11'
footprint = 'circle' # circle, rectangle or camera
ephMode = 'nbody' # nbody or 2body
ephType = 'basic' # basic (subset of ephemeride values - no helio, topo, etc.) or full (all that oorb calculates)
tStep = 1.0
roughTol = 10.0
obsFile = '_'.join([opsimrun, orbitfile.replace('.des', '').replace('.s3m', ''), 'obs.txt'])
obs = mo.DirectObs(footprint, rFov=rFov, tstep=tStep, roughTol=roughTol,
                   obsCode=obscode, obsTimeScale='TAI', ephMode=ephMode, ephType=ephType,
                   obsTimeCol=colmap['mjd'], obsRA=colmap['ra'], obsDec=colmap['dec'],
                   seeingCol=colmap['seeingGeom'], visitExpTimeCol=colmap['exptime'],
                   outfileName=obsFile, obsMetadata='%s test (night<100)' % (opsimrun))

[12]:

# Read the filter curves (to calculate the colors for each object)
filterlist = np.unique(simdata[colmap['filter']])
obs.readFilters(filterlist=filterlist)
# Calculate all colors ahead of time - this is just to create a dictionary early.
sednames = np.unique(orbits.orbits['sed_filename'])
for sedname in sednames:
    obs.calcColors(sedname)
# Want to see the lsst colors?  (these are all x-V band colors)
obs.colors

[12]:

{'C.dat': {'g': 0.28041793225569478,
  'i': -0.2926196728132453,
  'r': -0.177293845688272,
  'u': 1.5442111906705982,
  'y': -0.3028771664514025,
  'z': -0.29808899160965296},
 'S.dat': {'g': 0.36850606460945201,
  'i': -0.45720772665728049,
  'r': -0.26270911487248227,
  'u': 1.8626393418302385,
  'y': -0.4073469917610737,
  'z': -0.39986031775197262}}

[13]:

# Generate the actual observations and write to output file.
print('Will write output to %s' % obs.outfileName)
obs.run(orbits, simdata)

Will write output to baseline2018b_neos_100_obs.txt
Generating preliminary ephemerides on a grid of 1.000000 day timesteps.

[14]:

# This is what the output file looks like:
# it contains some metadata about how it was generated, and then the observations.
!head -16 $obsFile

# baseline2018b test (night<100)
# baseline2018b_neos_100_obs.txt
# ephemeris generation via PyOrbEphemerides
# planetary ephemeris file /Users/lynnej/lsstRepos/oorb/data/de430.dat
# obscode I11
# direct obs header metadata
# observation generation via DirectObs
# ephMode nbody prelimEphMode 2body
# rough tolerance for preliminary match 10.000000
# time step for preliminary match 1.000000
# pointing footprint circle
# rfov 1.750000
# obsRA fieldRA obsDec fieldDec obsRotSkyPos rotSkyPos obsDeg True
# obsMJD observationStartMJD obsTimeScale TAI seeing seeingFwhmGeom expTime visitExposureTime
objId time ra dec dradt ddecdt phase solarelon helio_dist geo_dist magV trueAnomaly velocity observationStartMJD night fieldRA fieldDec rotSkyPos filter visitExposureTime seeingFwhmEff seeingFwhmGeom fiveSigmaDepth solarElong magFilter dmagColor dmagTrail dmagDetect
S1000001 59854.0760648 335.661572918 -4.66872150058 -0.127477289493 -0.0991730463403 12.9028682217 146.83925488 2.45213370361 1.55215427431 21.6707534809 -99.0 12.9032493442 59854.0760648 2 335.151362 -5.689996 228.530394116 y 30.0 0.829128431145 0.733543570401 22.0670183628 145.792553804 21.2634064892 -0.407346991761 1.46493584452 2.8526839965

[15]:

# Using the "DirectObs" class.
rFov = 1.75  # in degrees
obscode = 'I11'
footprint = 'circle' # circle, rectangle or camera
ephMode = 'nbody' # nbody or 2body
ephType = 'full' # basic (subset of ephemeride values - no helio, topo, etc.) or full (all that oorb calculates)
tStep = 1.0
roughTol = 10.0
obsFile = '_'.join([opsimrun, orbitfile.replace('.des', '').replace('.s3m', ''), 'full_obs.txt'])
obs = mo.DirectObs(footprint, rFov=rFov, tstep=tStep, roughTol=roughTol,
                   obsCode=obscode, obsTimeScale='TAI', ephMode=ephMode, ephType=ephType,
                   obsTimeCol=colmap['mjd'], obsRA=colmap['ra'], obsDec=colmap['dec'],
                   seeingCol=colmap['seeingGeom'], visitExpTimeCol=colmap['exptime'],
                   outfileName=obsFile, obsMetadata='%s test (night<100)' % (opsimrun))

[16]:

# Read the filter curves (to calculate the colors for each object)
filterlist = np.unique(simdata[colmap['filter']])
obs.readFilters(filterlist=filterlist)
# Calculate all colors ahead of time - this is just to create a dictionary early.
sednames = np.unique(orbits.orbits['sed_filename'])
for sedname in sednames:
    obs.calcColors(sedname)

[17]:

# Generate the actual observations and write to output file.
print('Will write output to %s' % obs.outfileName)
obs.run(orbits, simdata)

Will write output to baseline2018b_neos_100_full_obs.txt
Generating preliminary ephemerides on a grid of 1.000000 day timesteps.

[18]:

# This is what the output file looks like:
!head -16 $obsFile

# baseline2018b test (night<100)
# baseline2018b_neos_100_full_obs.txt
# ephemeris generation via PyOrbEphemerides
# planetary ephemeris file /Users/lynnej/lsstRepos/oorb/data/de430.dat
# obscode I11
# direct obs header metadata
# observation generation via DirectObs
# ephMode nbody prelimEphMode 2body
# rough tolerance for preliminary match 10.000000
# time step for preliminary match 1.000000
# pointing footprint circle
# rfov 1.750000
# obsRA fieldRA obsDec fieldDec obsRotSkyPos rotSkyPos obsDeg True
# obsMJD observationStartMJD obsTimeScale TAI seeing seeingFwhmGeom expTime visitExposureTime
objId time ra dec dradt ddecdt phase solarelon helio_dist geo_dist magV pa topo_lon topo_lat opp_topo_lon opp_topo_lat helio_lon helio_lat opp_helio_lon opp_helio_lat topo_obj_alt topo_solar_alt topo_lunar_alt lunar_phase lunar_dist helio_x helio_y helio_z helio_dx helio_dy helio_dz obs_helio_x obs_helio_y obs_helio_z trueAnom velocity observationStartMJD night fieldRA fieldDec rotSkyPos filter visitExposureTime seeingFwhmEff seeingFwhmGeom fiveSigmaDepth solarElong magFilter dmagColor dmagTrail dmagDetect
S1000001 59854.0760648 335.661572918 -4.66872150058 -0.127477289493 -0.0991730463403 12.9028682217 146.840912582 2.45213370361 1.55215427431 21.6707534809 232.118283673 335.741135515 5.08952875392 -32.810127053 5.09083443774 348.541969683 3.21849979771 -20.0092928851 3.21980548153 63.1861602278 -38.4665724493 36.1010079827 0.38739595815 70.6520986767 2.3994731612 -0.486348214868 0.137672302874 0.00218415563495 0.0103172709944 -0.000681661036841 0.989955765886 0.148855498807 -2.28131980833e-05 -99.0 12.9032493442 59854.0760648 2 335.151362 -5.689996 228.530394116 y 30.0 0.829128431145 0.733543570401 22.0670183628 145.792553804 21.2634064892 -0.407346991761 1.46493584452 2.8526839965

[19]:

# Or nicely read into a pandas dataframe..
d = pd.read_csv(obsFile, delim_whitespace=True, comment='#')
d[0:5]

[19]:
objId time ra dec dradt ddecdt phase solarelon helio_dist geo_dist ... filter visitExposureTime seeingFwhmEff seeingFwhmGeom fiveSigmaDepth solarElong magFilter dmagColor dmagTrail dmagDetect
0 S1000001 59854.076065 335.661573 -4.668722 -0.127477 -0.099173 12.902868 146.840913 2.452134 1.552154 ... y 30.0 0.829128 0.733544 22.067018 145.792554 21.263406 -0.407347 1.464936 2.852684
1 S1000001 59854.083738 335.660591 -4.669482 -0.127479 -0.099123 12.905885 146.832261 2.452134 1.552209 ... y 30.0 0.777921 0.691451 22.123200 147.792759 21.263484 -0.407347 1.496139 2.914742
2 S1000001 59856.017419 335.442640 -4.856711 -0.113656 -0.094872 13.637453 144.697815 2.452229 1.566566 ... z 30.0 1.015172 0.886471 22.712273 143.899147 21.290969 -0.399860 1.394707 2.712270
3 S1000001 59856.038565 335.440224 -4.858716 -0.113971 -0.094735 13.645485 144.674272 2.452230 1.566726 ... z 30.0 0.786676 0.698648 23.004635 143.878362 21.291190 -0.399860 1.520039 2.962138
4 S1000001 59861.079375 334.989973 -5.303635 -0.079843 -0.081813 15.435813 139.204684 2.452413 1.607632 ... z 30.0 0.699799 0.627234 22.641844 138.945878 21.347170 -0.399860 1.642281 3.202670

5 rows × 51 columns

[20]:

# Note that the apparent magnitude in any given visit will be generated by MAF, or can be calculated as:
d['appMag'] = d['magV'] + d['dmagColor'] + d['dmagDetect']
cols = ['objId', 'magV', 'dmagColor',  'magFilter', 'filter',
        'seeingFwhmGeom', 'dradt', 'ddecdt', 'dmagDetect',
        'appMag', 'fiveSigmaDepth']
tmp = d[cols][23:33]
tmp

[20]:
objId magV dmagColor magFilter filter seeingFwhmGeom dradt ddecdt dmagDetect appMag fiveSigmaDepth
23 S1000001 22.114624 -0.399860 21.714764 z 0.885394 0.108423 -0.002984 3.215100 24.929864 22.400191
24 S1000001 22.114850 -0.399860 21.714990 z 0.936405 0.108590 -0.002860 3.157091 24.872081 22.395591
25 S1000001 22.142264 -0.399860 21.742403 z 0.714136 0.120133 0.002469 3.449325 25.191728 22.631871
26 S1000001 22.142476 -0.399860 21.742615 z 0.628592 0.120305 0.002583 3.580050 25.322666 22.909896
27 S1000004 12.801211 -0.302877 12.498334 y 0.787645 0.339752 -0.097434 3.352014 15.850348 21.768465
28 S1000005 9.375908 -0.302877 9.073030 y 0.717491 0.159049 -0.250901 3.661206 12.734236 21.991805
29 S1000005 9.375825 -0.302877 9.072947 y 0.637436 0.158952 -0.250797 3.780992 12.853939 21.993364
30 S1000005 9.262223 1.544211 10.806434 u 0.626873 0.110650 -0.235304 3.779455 14.585889 23.613301
31 S1000005 9.247978 1.544211 10.792189 u 0.560244 0.104324 -0.232891 3.889496 14.681685 23.786732
32 S1000005 9.176317 -0.302877 8.873440 y 0.718928 0.071651 -0.218266 3.620191 12.493631 22.121915 
[21]:

orbits.orbits.query('objId in @tmp.objId')

[21]:
objId otype q e inc Omega argPeri tPeri H g epoch model sed_filename
1 S1000001 S3M 2.10974 0.07518 4.91571 209.40298 322.66447 56722.78779 20.57 0.15 57231.0 V150728 S.dat
4 S1000004 S3M 2.17676 0.19949 12.92422 162.14580 192.22312 56808.90252 10.56 0.15 57231.0 V150728 C.dat
5 S1000005 S3M 1.98262 0.16324 23.42404 210.02222 326.18514 57395.29592 8.17 0.15 57231.0 V150728 C.dat 
[ ]: