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

In [1]:

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

In [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.¶

In [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))

In [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']

In [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'.

In [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]

Out[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.¶

In [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)

In [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']

In [9]:

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

In [10]:

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

Out[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).

In [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))

In [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

Out[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}}

In [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.

In [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

In [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))

In [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)

In [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.

In [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

In [19]:

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

Out[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

In [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

Out[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

In [21]:

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

Out[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