Sphinx logo

obstools – miscellaneous tools for observation

The obstools module stores tools for oberving (pre- and post-) as well as functioning as a module for various corrections and simple calculations that don’t have a better place to live.

Most implementations are for optical astronomy, as that is what the primary author does.

Note that throughout this module it is assumed that UTC as used in datetime is the same as UT1 (the UT used for calculations). This is only an issue if leapseconds stop being updated or if something such as daylight savings time calculations in datetime are important to be correct to less than DUT1 - otherwise, the “UTC” in datetime can simply by replaced with UT1 values when precision better than DUT1 is needed.

Also note that some of these functions require the dateutil package (it is installed by default with matplotlib)

See also

astropysics.coords.earth_rotation_angle() astropysics.coords.greenwich_sidereal_time()

Classes and Inheritance Structure

Inheritance diagram of astropysics.obstools

Module API

class astropysics.obstools.CalzettiExtinction(A0=1)

Bases: astropysics.obstools.Extinction

This is the average extinction law derived in Calzetti et al. 1994: x=1/lambda in mu^-1 Q(x)=-2.156+1.509*x-0.198*x**2+0.011*x**3

f(lamb)
class astropysics.obstools.CardelliExtinction(EBmV=1, Rv=3.1)

Bases: astropysics.obstools._EBmVExtinction

Milky Way Extinction law from Cardelli et al. 1989

f(lamb)
class astropysics.obstools.Extinction(f=None, A0=1)

Bases: object

This is the base class for extinction-law objects. Extinction laws can be passed in as a function to the initializer, or subclasses should override the function f with the preferred law as f(self,lambda), and their initializers should set the zero point

Note that functions are interpreted as magnitude extinction laws ... if optical depth is desired, f should return (-2.5/log(10))*tau(lambda)

A0 is the normalization factor that gets multiplied into the reddening law.

Alambda(band)

determines the extinction for this extinction law in a given band or bands

band can be a wavelength or a string specifying a band

computeA0FromFluxRatio(measured, expected, lambda1=None, lambda2=None, filterfunc=None)

This derives the normalization of the Extinction function from provided ratios for theoretically expected fluxes. If multiple measurements are provided, the mean will be used

measured is the measured line ratio (possibly an array), while expected is either the expected line ratios, or a string specifying the appropriate balmer flux ratio as “Hab”,”Hde”,etc. (for Halpha/Hbeta or Hdelta/Hepsilon) to assume case B recombination fluxes (see e.g. Osterbrock 2007).

lambda1 and lambda2 are the wavelengths for the ratios F1/F2, or None if a string is provided

filterfunc is a function to be applied as the last step - it can either be used to contract an array to (e.g. np.mean), or filter out invalid values (e.g. lambda x:x[np.isfinite(x)]). Default does nothing

returns A0,standard deviation of measurements

correctColor(colors, bands)

Uses the supplied extinction law to correct a color (or array of colors) where the color is in the specified bands

bands is a length-2 sequence with either a band name or a wavelength for the band, or of the form ‘bandname1-bandname2’ or ‘E(band1-band2)’

correctFlux(flux, lamb)
correctPhotometry(mags, band)

Uses the extinction law to correct a magnitude (or array of magnitudes)

bands is either a string specifying the band, or a wavelength to use

correctSpectrum(spec, newspec=True)

Uses the supplied extinction law to correct a spectrum for extinction.

if newspec is True, a copy of the supplied spectrum will have the extinction correction applied

returns the corrected spectrum

f(lamb)
class astropysics.obstools.FMExtinction(C1, C2, C3, C4, x0, gamma, EBmV=1, Rv=3.1)

Bases: astropysics.obstools._EBmVExtinction

Base class for Extinction classes that use the form from Fitzpatrick & Massa 90

f(lamb)
class astropysics.obstools.LMCExtinction(EBmV=0.3, Rv=3.41)

Bases: astropysics.obstools.FMExtinction

LMC Extinction law from Gordon et al. 2003 LMC Average Sample

class astropysics.obstools.SMCExtinction(EBmV=0.2, Rv=2.74)

Bases: astropysics.obstools.FMExtinction

SMC Extinction law from Gordon et al. 2003 SMC Bar Sample

class astropysics.obstools.Site(lat, long, alt=0, tz=None, name=None)

Bases: object

This class represents a location on Earth from which skies are observable.

lat/long are coords.AngularCoordinate objects, altitude in meters. latitude here is the geographic/geodetic latitude.

generate a site by specifying latitude and longitude (as coords.AngularCoordinate objects or initializers for one), optionally providing altitude (in meters), time zone (either as a timezone name provided by the system or as an offset from UTC), and/or a site name.

altitude

Altitude of the site in meters

apparentCoordinates(coords, datetime=None, precess=True, refraction=True)

computes the positions in horizontal coordinates of an object with the provided fixed coordinates at the requested time(s).

datetime can be a sequence of datetime objects or similar representations, or it can be a sequnce of JD’s. If None, it uses the current time or the value of the currentobsjd property.

If precess is True, the position will be precessed to the epoch of the observation (almost always the right thing to do)

If refraction is True, an added correction to the altitude due to atmospheric refraction at STP (formula from Meeus ch 16) is included. If refraction is a (non-0) float, it will be taken as the temperature at which to perform the refraction calculation. If it evaluates to False, no refraction correction is performed

returns a sequence of astropysics.coords.HorizontalCoordinates objects

currentobsjd

Date and time to use for computing time-dependent values. If set to None, the jd at the instant of calling will be used. It can also be set as datetime objects or (yr,mon,day,hr,min,sec) tuples.

equatorialToHorizontal(eqpos, lsts, epoch=None)

Generates a list of horizontal positions (or just one) from a provided equatorial position and local siderial time (or sequence of LSTs) in decimal hours and a latitude in degrees. If epoch is not None, it will be used to set the epoch in the equatorial system.

Note that this generally should not be used to compute the observed horizontal coordinates directly - apparentCoordinates() includes all the corrections and formatting. This method is purely for the coordinate conversion.

geocentriclat
latitude

Geographic/Geodetic Latitude of the site as an AngularCoordinate object

localSiderialTime(*args, **kwargs)

Compute the local siderial time given an input civil time. The various input forms are used to determine the interpretation of the civil time:

  • localSiderialTime()

    current local siderial time for this Site or uses the value of the currentobsjd property.

  • localSiderialTime(JD)

    input argument is julian date UT1

  • localSdierialTime(datetime.date)

    compute the local siderial time for midnight on the given date

  • localSiderialTime(datetime.datetime)

    the datetime object specifies the local time. If it has tzinfo, the object’s time zone will be used, otherwise the Site's

  • localSiderialTime(time,year,month,day)

    input arguments determine local time - time is in hours

  • localSiderialTime(year,month,day,hr,min,sec)

    local time - hours and minutes will be interpreted as integers

keywords

  • apparent

    if True (default) the returned time will be local apparent sidereal time (i.e. nutation terms included), otherwise it will be local mean sidereal time

  • returntype

    a string that determines the form of the returned LST as described below

returns the local siderial time in a format that depends on the returntype keyword. It can be:

  • None/’hours’ (default)

    LST in decimal hours

  • ‘string’

    LST as a hh:mm:ss.s

  • ‘datetime’

    a datetime.time object

localTime(lsts, date=None, apparent=True, returntype=None, utc=False)

Computes the local civil time given a particular local siderial time.

lsts are the input local times, and may be either a scalar or array, and must be in decimal hours. Althernatively, it can be a datetime.datetime object, in which case the date will be inferred from this object and the date argument will be ignored.

date should be a datetime.date object or a (year,month,day) tuple. If None, the current date will be assumed as inferred from Site.currentobsjd

if lsts is a datetime.datetime object, the object will be interpreted as the local siderial time with the corresponding date (any timezone info will be ignored), and the date argument will be ignored.

if apparent is True, the inputs are assumed to be in local apparent siderial time (i.e. nutation terms included), otherwise local mean siderial time.

returns the local time in a format that depends on the returntype keyword. It can be:

  • None/’hours’ (default)

    local time in decimal hours

  • ‘string’

    local time as a hh:mm:ss.s

  • ‘datetime’

    a datetime.time object with the appropriate tzinfo

If utc is True, the time will be converted to UTC before being returned.

longitude

Longitude of the site as an AngularCoordinate object

nextRiseSetTransit(eqpos, dtime=None, alt=-0.5667)

Get ‘next’ rise,set,transit for an equatorial position at a given datetime (see details).

If on-sky at the specified datetime, returns the rise / set times bracketing current transit, i.e. returns rise, set, transit such that rise < now < set; so rise and possibly transit will be in the past at specified time.

If the target is off-sky currently, then rise/transit/set will all be in the future, so now < rise < transit < set.

If circumpolar (always visible) rise/set are None, and next transit is returned. If never visible, rise / set / transit are all None.

NB The only dates checked for visibility are ‘today’ and ‘tomorrow’ (relative to given datetime). If the target is an edge case that’s below the horizon all day today, but is visible in 6 months, you will still get a None result.

See also: riseSetTransit

args alt` determines the altitude to be considered as risen or set in degrees. Default is for approximate rise/set including refraction.

dtime determines the datetime at which to do the computation. See calendar_to_jd() for acceptable formats. If None, the current datetime will be assumed as inferred from Site.currentobsjd

returns (rise,set,transit) as :class:datetime.datetime objects (UTC). If the object is circumpolar, rise and set are both None. If it is never visible, rise, set, and transit are all None.

nightPlot(coords, date=None, plottype='altam', onlynight=False, clf=True, utc=False, colors=None, plotkwargs=None)

Generates plots of important observability quantities for the provided coordinate objects for a single night.

coords should be a astropysics.coords.LatLongCoordinates, a sequence of such objects, or a dictionary mapping the name of an object to the object itself. These names will be used to label the plot.

plottype can be one of:

  • ‘altam’ : a plot of time vs. altitude with a secondary axis for sec(z)
  • ‘am’ : a plot of time vs. airmass/sec(z)
  • ‘altaz’: a plot of azimuth vs. altitude
  • ‘sky’: polar projection of the path on the sky

If onlynight is True, the plot will be shrunk to only show the times between sunset and sunrise. Otherwise, an entire day/night will be plotted.

If clf is True, the figure is cleared before the observing plot is made.

If utc is True, the times will be in UTC instead of local time.

colors should be a sequence of matplotlib color specifiers, or None to use the default color cycle.

plotkwargs will be provided as a keyword dictionary to matplotlib.pyplot.plot(), unless it is None

nightTable(coord, date=None, strtablename=None, localtime=True, hrrange=(18, 6, 25))

tabulates altitude, azimuth, and airmass values for the provided fixed position on a particular date, specified either as a datetime.date object, a (yr,mon,day) tuple, or a julain date (will be rounded to nearest)

if date is None, the current date is used

If localtime is True, the hours output (and input) will be in local time for this Site. Otherwise, it is UTC.

hrrange determines the size of the table - it should be a 3-tuple (starthr,endhr,n) where starthr is on the day specified in date and endhr is date + 1 day

if strtablename is True, a string is returned with a printable table of observing data. If strtable is a string, it will be used as the title for the table. Otherwise, a record array is returned with the hour(UTC), alt, az, and airmass

onSky(eqpos, dtime=None, alt=-0.5667)

Checks if a target is on-sky at specified datetime.

NB on-sky means above minimum elevation as specifed by ‘alt’ parameter.

args eqpos Ra,Dec in equatorial co-ordinates.

dtime determines the datetime at which to do the computation. See calendar_to_jd() for acceptable formats. If None, the current datetime will be assumed as inferred from Site.currentobsjd

alt determines the altitude to be considered as risen or set in degrees. Default is for approximate rise/set including refraction.

returns
Boolean.
riseSetTransit(eqpos, date=None, alt=-0.5667, timeobj=False, utc=False)

Computes the rise, set, and transit times of a provided equatorial position in local time.

alt determines the altitude to be considered as risen or set in degrees. Default is for approximate rise/set including refraction.

date determines the date at which to do the computation. See calendar_to_jd() for acceptable formats. Note that this date is the date of the transit, so the rise and set times may be on the previous/following day.

returns (rise,set,transit) as :class:datetime.datetime objects if timeobj is True or if it is False, they are in decimal hours. If the object is circumpolar, rise and set are both None. If it is never visible, rise,set, and transit are all None

yearPlot(coords, startdate=None, n=13, months=12, utc=False, clf=True, colors=None)

Plots the location transit and rise/set times of object(s) over ~year timescales.

coords should be a astropysics.coords.LatLongCoordinates, a sequence of such objects, or a dictionary mapping the name of an object to the object itself. These names will be used to label the plot.

startdate is a datetime.date object or a date tuple (year,month,day) that is used as the start of the plot.

n specifies the number of points to include of the objects

months is the number of months to draw the plot for.

If utc is True, the times will be in UTC instead of local time.

If clf is True, the figure is cleared before the observing plot is made.

colors should be a sequence of matplotlib color specifiers, or None to use the default color cycle.

astropysics.obstools.calendar_to_jd(caltime, tz=None, gregorian=True, mjd=False)

Convert a calendar date and time to julian date.

Parameters:
  • caltime

    The date and time to compute the JD. Can be in one of these forms:

    • A sequence of floats in the order (yr,month,day,[hr,min,sec]).
    • A sequence in the order (yr,month,day,[hr,min,sec]) where at least
      one of the elements is a sequence (a sequence will be returned).
    • A datetime.datetime or datetime.date object
    • A sequence of datetime.datetime or datetime.date objects (a sequence will be returned).
    • None : returns the JD at the moment the function is called.

    If the time is unspecified, it is taken to be noon (i.e. Julian Date = Julian Day Number)

  • tz

    Sets the time zone to assume for the inputs for conversion to UTC. Can be any of the following:

    • None
      No time zone conversion will occur unless caltime is given as datetime.datetime or datetime.date objects with tzinfo, in which case they will be converted to UTC using their own tzinfo.
    • a string
      Specifies a timezone name (resolved into a timezone using the dateutil.tz.gettz() function).
    • a scalar
      The hour offset of the timezone.
    • a datetime.tzinfo object,
      This object will be used for timezone information.
  • gregorian (bool or None) – If True, the input will be interpreted as in the Gregorian calendar. Otherwise, it will be Julian. If None, it will be assumed to switch over on October 4/15, 1582.
  • mjd (bool) – If True, a modified julian date is returned instead of the standard julian date.
Returns:

JD as a float, or a sequence of JDs if sequences were input.

Examples

>>> import datetime,dateutil
>>> calendar_to_jd((2010,1,1))
2455198.0
>>> calendar_to_jd(datetime.datetime(2000,12,21,3,0,0))
2451899.625
>>> calendar_to_jd([2004,3,(5,6)])
array([ 2453070.,  2453071.])
>>> dates = [datetime.datetime(2004,3,5),datetime.datetime(2004,3,9)]
>>> calendar_to_jd(dates)
array([ 2453069.5,  2453073.5])
>>> tz = dateutil.tz.tzoffset('2',3*3600)
>>> calendar_to_jd((2010,1,1),tz)
2455197.875
astropysics.obstools.delta_AT(jdutc, usett=False)

Computes the difference between International Atomic Time (TAI) and UTC, known as delta(AT).

Note that this is not valid before UTC (Jan 1,1960) beganand it is not correct for future dates, as leap seconds are not predictable. Hence, warnings are issued if before UTC or >5 years from the date of the algorithm.

Implementation adapted from the matching SOFA algorithm (dat.c).

Parameters:
  • utc (float) – UTC time as a Julian Date (use calendar_to_jd() for calendar form inputs.)
  • usett (bool) – If True, the return value will be the difference between UTC and Terrestrial Time (TT) instead of TAI (TT - TAI = 32.184 s).
Returns:

TAI - UTC in seconds as a float (or TT - UTC if usett is True)

astropysics.obstools.epoch_to_jd(epoch, julian=True, mjd=False)

Converts a Julian or Besselian Epoch to a Julian Day.

Parameters:
  • epoch (string, scalar, or array-like) – The epoch as a decimal year.
  • julian (bool) – If True, a Julian Epoch will be used (the year is exactly 365.25 days long). Otherwise, the epoch will be Besselian (assuming a tropical year of 365.242198781 days). If epoch is a string and starts with ‘B’ or ‘J’, this parameter will be ignored, and the ‘B’ or ‘J’ specifies the epoch type.
  • mjd (bool) – If True, a modified julian date is returned instead of the standard julian date.
Returns:

The Julian Day as a float or array (if epoch is array-like)

Reference:

http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epj.html

astropysics.obstools.extinction_correction(lineflux, linewl, EBmV, Rv=3.1, exttype='MW')

Extinction correct a la Cardelli et al 89 from the supplied line data and a given E(B-V) along the sightline

inputs may be numpy arrays

linewl is in angstroms, lineflux in erg s^-1 cm^-2

if the input lineflux is None (or NaN but NOT False or 0) , Alambda is returned instead

astropysics.obstools.extinction_from_flux_ratio(frobs, frexpect, outlambda=None, Rv=3.1, tol=0.0001)

frobs is the observed flux ratio f1/f2

frexpect is either a string code specifying a hydrogen transition (e.g. ‘Hab’ is Halpha/Hbeta, ‘Hde’ is Hdelta/Hepsilon, from Ostriker 2E) or a tuple of the form (expected f1/f2,lambda1,lambda2) wl in angstroms

outputs E(B-V) to a tolerance specified by tol if outlambda is 0/False/None, otherwise outputs Alambda (tol still determines E(B-V)) (outlambda can be UBVRI or ugriz as from B&M)

frobs can be an array, but other values cannot

astropysics.obstools.get_SFD_dust(long, lat, dustmap='ebv', interpolate=True)

Gets map values from Schlegel, Finkbeiner, and Davis 1998 extinction maps.

dustmap can either be a filename (if ‘%s’ appears in the string, it will be replaced with ‘ngp’ or ‘sgp’), or one of:

  • ‘i100’

    100-micron map in MJy/Sr

  • ‘x’

    X-map, temperature-correction factor

  • ‘t’

    Temperature map in degrees Kelvin for n=2 emissivity

  • ‘ebv’

    E(B-V) in magnitudes

  • ‘mask’

    Mask values

For these forms, the files are assumed to lie in the current directory.

Input coordinates are in degrees of galactic latiude and logitude - they can be scalars or arrays.

if interpolate is an integer, it can be used to specify the order of the interpolating polynomial

astropysics.obstools.get_dust_radec(ra, dec, dustmap, interpolate=True)
astropysics.obstools.jd_to_calendar(jd, rounding=1000000, output='datetime', gregorian=None, mjd=False)

Converts a julian date to a calendar date and time.

Parameters:
  • jd (scalar, array-like, or None) – The Julian Date at which to compute the calendar date/time, a sequence of JDs, or None for the current date/time at the moment the function is called.
  • rounding (scalar) – If non-0, Performs a fix for floating-point errors. It specifies the number of milliseconds by which to round the result to the nearest second. If 1000000 (one second), no milliseconds are recorded. If larger, a ValueError is raised.
  • output

    Determines the format of the returned object and can be:

    • ‘datetime’
      A list of datetime.datetime objects in UTC will be returned. If the input is a scalar, a single object will be returned.
    • ‘array’
      A Nx7 array will be returned of the form [(year,month,day,hr,min,sec,msec),...] unless the input was a scalar, in which case it will be a length-7 array.
    • ‘fracarray’
      An Nx3 array (year,month,day) where day includes the decimal portion.
  • gregorian (bool or None) – If True, the output will be in the Gregorian calendar. Otherwise, it will be Julian. If None, it will be assumed to switch over on October 4/15 1582.
  • mjd (bool) – If True, the input is interpreted as a modified julian date instead of a standard julian date.
Returns:

The calendar date and time in a format determined by the output parameter (see above).

Raises ValueError:
 

If rounding is larger than one second, or output is invalid.

Examples

>>> jd_to_calendar(2451545)
datetime.datetime(2000, 1, 1, 12, 0, tzinfo=tzutc())
>>> jd_to_calendar(2305812.5)
datetime.datetime(1600, 12, 31, 0, 0, tzinfo=tzutc())
>>> jd_to_calendar([2415020.5,2305447.5],output='array')
array([[1900,    1,    1,    0,    0,    0,    0],
       [1600,    1,    1,    0,    0,    0,    0]])
>>> jd_to_calendar(0.0,output='fracarray')
array([[ -4.71200000e+03,   1.00000000e+00,   1.50000000e+00]])
astropysics.obstools.jd_to_epoch(jd, julian=True, asstring=False, mjd=False)

Converts a Julian Date to a Julian or Besselian Epoch expressed in decimal years.

Parameters:
  • jd (scalar, array-like, or None) – Julian Date for computing the epoch, or None for current epoch.
  • julian (bool) – If True, a Julian Epoch will be used (the year is exactly 365.25 days long). Otherwise, the epoch will be Besselian (assuming a tropical year of 365.242198781 days).
  • asstring (bool) – If True, a string of the form ‘J2000.0’ will be returned. If it is an integer, the number sets the number of significant figures in the output string Otherwise, scalars are returned (an int if a whole year, float if not).
  • mjd (bool) – If True, a modified julian date is returned instead of the standard julian date.
Returns:

The epoch as a string (or list of strings if jd was array-like) if asstring is True. If not, an int if a whole year, or a float (or array if jd was array-like).

Reference:

http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epj.html

astropysics.obstools.mjdoffset = 2400000.5

Offset between Julian Date and Modified Julian Date - e.g. mjd = jd - mjdoffset