Source code for

# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-

I/O utility functions for files in desimodel.
import os
import re
import warnings
from datetime import datetime, timezone

import yaml
import gzip
import numpy as np
from import fits
from astropy.table import Table, Column

from desiutil.log import get_logger

log = get_logger()

_thru = dict()

# ADM raise a custom exception when an environment variable is missing.
[docs]class MissingEnvVar(Exception): pass
[docs]def load_throughput(channel): """Returns specter Throughput object for the given channel 'b', 'r', or 'z'. Parameters ---------- channel : {'b', 'r', 'z'} Spectrograph channel. Returns ------- Throughput A specter throughput object. """ import specter.throughput channel = channel.lower() global _thru if channel not in _thru: thrufile = findfile("throughput/thru-{0}.fits".format(channel)) _thru[channel] = specter.throughput.load_throughput(thrufile) return _thru[channel]
# # # _psf = dict()
[docs]def load_psf(channel): """Returns specter PSF object for the given channel 'b', 'r', or 'z'. Parameters ---------- channel : {'b', 'r', 'z'} Spectrograph channel. Returns ------- PSF A specter PSF object. """ import specter.psf channel = channel.lower() global _psf if channel not in _psf: psffile = findfile("specpsf/psf-{0}.fits".format(channel)) _psf[channel] = specter.psf.load_psf(psffile) return _psf[channel]
# # # _params = None
[docs]def load_desiparams(): """Returns DESI parameter dictionary loaded from ``$DESIMODEL/data/desi.yaml``. Returns ------- :class:`dict` The parameters read from the YAML file. """ global _params if _params is None: desiparamsfile = findfile("desi.yaml") with open(desiparamsfile) as par: _params = yaml.safe_load(par) # - for temporary backwards compability after 'exptime' -> 'exptime_dark' if ("exptime" not in _params) and ("exptime_dark" in _params): _params["exptime"] = _params["exptime_dark"] # - Augment params with wavelength coverage from specpsf files # - wavemin/max = min/max wavelength covered by *any* fiber on the CCD # - wavemin/max_all = min/max wavelength covered by *all* fibers for channel in ["b", "r", "z"]: hdr = fits.getheader(findfile("specpsf/psf-{}.fits".format(channel)), 0) _params["ccd"][channel]["wavemin"] = hdr["WAVEMIN"] _params["ccd"][channel]["wavemax"] = hdr["WAVEMAX"] _params["ccd"][channel]["wavemin_all"] = hdr["WMIN_ALL"] _params["ccd"][channel]["wavemax_all"] = hdr["WMAX_ALL"] return _params
# # # # Added and still needs to be committed and pushed to desihub _gfa = None
[docs]def load_gfa(): """Returns GFA table from ``$DESIMODEL/data/focalplane/gfa.ecsv``. Returns ------- :class:`~astropy.table.Table` The data from the ECSV file. """ global _gfa if _gfa is None: gfaFile = findfile("focalplane/gfa.ecsv") _gfa =, format="ascii.ecsv") return _gfa
# # # _deviceloc = None
[docs]def load_deviceloc(): """Returns a table from ``$DESIMODEL/data/focalplane/fiberpos-all.fits``. Returns ------- :class:`~astropy.table.Table` The data from the FITS file, with columns converted to uppercase. """ global _deviceloc if _deviceloc is None: fiberposfile = findfile("focalplane/fiberpos-all.fits") _deviceloc = # - Convert to upper case if needed # - Make copy of colnames b/c they are updated during iteration for col in list(_deviceloc.colnames): if col.islower(): _deviceloc.rename_column(col, col.upper()) return _deviceloc
_fiberpos = None
[docs]def load_fiberpos(): """Returns fiberpos table from ``$DESIMODEL/data/focalplane/fiberpos.fits``. Returns ------- :class:`~astropy.table.Table` The data from the FITS file, sorted by ``FIBER``. """ global _fiberpos if _fiberpos is None: fiberposfile = findfile("focalplane/fiberpos.fits") _fiberpos = _fiberpos.sort("FIBER") # - Convert to upper case if needed # - Make copy of colnames b/c they are updated during iteration for col in list(_fiberpos.colnames): if col.islower(): _fiberpos.rename_column(col, col.upper()) # - Temporary backwards compatibility for renamed columns if "POSITIONER" in _fiberpos.colnames: warnings.warn( "old fiberpos.fits with POSITIONER column instead of LOCATION; please update your $DESIMODEL checkout", DeprecationWarning, ) _fiberpos["LOCATION"] = _fiberpos["POSITIONER"] else: _fiberpos["POSITIONER"] = _fiberpos["LOCATION"] if "SPECTROGRAPH" in _fiberpos.colnames: warnings.warn( "old fiberpos.fits with SPECTROGRAPH column instead of SPECTRO; please update your $DESIMODEL checkout", DeprecationWarning, ) _fiberpos["SPECTRO"] = _fiberpos["SPECTROGRAPH"] else: _fiberpos["SPECTROGRAPH"] = _fiberpos["SPECTRO"] return _fiberpos
# # # _tiles = dict()
[docs]def load_tiles(onlydesi=True, extra=False, tilesfile=None, cache=True, programs=None, surveyops=True): """Return DESI tiles structure from ``$DESI_SURVEYOPS/trunk/ops/tiles-main.ecsv``. Parameters ---------- onlydesi : :class:`bool`, optional If ``True``, trim to just the tiles in the DESI footprint. extra : :class:`bool`, optional If ``True``, include extra layers with ``PROGRAM='EXTRA'``. tilesfile : :class:`str`, optional Name of tiles file to load; or None for default. See Notes for details. cache : :class:`bool`, optional If ``False``, force reload of data from tiles file, instead of using cached values. programs : :class:`list` or `str`, optional Pass a list of program names to restrict to only those programs, e.g. ["DARK", "BACKUP"]. surveyops : :class:`bool` If ``True`` then find the relevant path for the $DESI_SURVEYOPS directory rather than the $DESIMODEL directory. Returns ------- :class:`` The data table portion of the FITS file. Raises ------ :exc:`FileNotFoundError` If the value of `tilesfile` does not exist. Notes ----- Keyword-based environment variable expansion is performed on the `tilesfile` value, so *e.g.*:: tiles = load_tiles(tilesfile='{HOME}/my-tiles.fits') will be expanded with the value of :envvar:`HOME`. If the parameter `tilesfile` is set, this function uses the following search method: 0. Paths corresponding to both $DESI_SURVEYOPS/trunk/ops and $DESI_SURVEYOPS/ops are always both checked, to cover different svn checkout approaches. 1. If the value includes an explicit path, even ``./``, use that file. 2. If the value does *not* include an explicit path, *and* the file name is identical to a file in ``$DESI_SURVEYOPS/trunk/ops/``, use the file in ``$DESI_SURVEYOPS/trunk/ops/`` and issue a warning. 3. If no matching file can be found at all, raise an exception. """ global _tiles if tilesfile is None: # Use the default if surveyops: tilesfile = findfile("tiles-main.ecsv", surveyops=surveyops) else: tilesfile = findfile("footprint/desi-tiles.fits", surveyops=surveyops) else: # If full path isn't included, check local vs $DESI_SURVEYOPS/ops tilepath, filename = os.path.split(tilesfile) if tilepath == "": have_local = os.path.isfile(tilesfile) checkfile = findfile(tilesfile, surveyops=surveyops) have_dmdata = os.path.isfile(checkfile) if have_dmdata: if have_local: msg = ( "$DESI_SURVEYOPS/(trunk)/ops/{0} is shadowed by a local" + " file. Choosing $DESI_SURVEYOPS file." + ' Use tilesfile="./{0}" if you want the local copy' + " instead." ).format(tilesfile) warnings.warn(msg) tilesfile = checkfile if not (have_local or have_dmdata): msg = ( 'File "{}" does not exist locally or in ' + "$DESI_SURVEYOPS/(trunk)/ops/!" ).format(tilesfile) raise FileNotFoundError(msg) # - standarize path location tilesfile = os.path.abspath(tilesfile.format(**os.environ)) log.debug("Loading tiles from %s", tilesfile) # ADM allow reading from either .fits or .ecsv files. # ADM guard against the possibility that the file is zipped. isfits = ".fits" in os.path.basename(tilesfile) if cache and tilesfile in _tiles: tiledata = _tiles[tilesfile] else: if isfits: with, memmap=False) as hdulist: tiledata = hdulist[1].data # # Temporary workaround for problem identified in # # if any([c.bzero is not None for c in tiledata.columns]): foo = [_tiles[k].dtype for k in tiledata.dtype.names] # - Check for out-of-date tiles file if np.issubdtype(tiledata["OBSCONDITIONS"].dtype, np.unsignedinteger): warnings.warn( "Old desi-tiles.fits with uint16 OBSCONDITIONS; please update your $DESIMODEL checkout.", DeprecationWarning, ) else: tiledata = # - load cache for next time if cache: _tiles[tilesfile] = tiledata # - Filter to only the DESI footprint if requested subset = np.ones(len(tiledata), dtype=bool) if onlydesi: subset &= tiledata["IN_DESI"] > 0 # - Filter out PROGRAM=EXTRA tiles if requested if not extra: subset &= ~np.char.startswith(tiledata["PROGRAM"], "EXTRA") # ADM filter to program names if requested. if programs is not None: # ADM guard against a single string being passed. programs = np.atleast_1d(programs) isprog = np.zeros(len(tiledata), dtype=bool) for program in programs: isprog |= tiledata["PROGRAM"] == program subset &= isprog if np.all(subset): return tiledata else: return tiledata[subset]
_platescale = None
[docs]def load_platescale(): """Loads platescale.txt. Returns ------- :class:`~numpy.recarray` The data table read from the file. Notes ----- The returned object has these columns: radius Radius from center of focal plane [mm]. theta Radial angle that has a centroid at this radius [deg]. radial_platescale Meridional (radial) plate scale [um/arcsec]. az_platescale: Sagittal (azimuthal) plate scale [um/arcsec]. arclength: Unknown description. """ global _platescale if _platescale is not None: return _platescale infile = findfile("focalplane/platescale.txt") columns = [ ("radius", "f8"), ("theta", "f8"), ("radial_platescale", "f8"), ("az_platescale", "f8"), ("arclength", "f8"), ] try: _platescale = np.loadtxt(infile, usecols=[0, 1, 6, 7, 8], dtype=columns) except (IndexError,ValueError): # - no "arclength" column in this version of desimodel/data # - Get info from separate rzs file instead _platescale = np.loadtxt(infile, usecols=[0, 1, 6, 7, 7], dtype=columns) rzs ="focalplane/rzsn.txt"), format="ascii") from scipy.interpolate import interp1d from numpy.lib.recfunctions import append_fields arclength = interp1d(rzs["R"], rzs["S"], kind="quadratic") _platescale["arclength"] = arclength(_platescale["radius"]) return _platescale
_focalplane = None def ensure_focalplane_loaded(): global _focalplane if _focalplane is not None: return # First call, parse all data files. fpdir = os.path.join(datadir(), "focalplane") fppat = re.compile(r"^desi-focalplane_(.*)\.ecsv$") stpat = re.compile(r"^desi-state_(.*)\.ecsv$") expat = re.compile(r"^desi-exclusion_(.*)\.(?:yaml|json).*$") fpraw = dict() msg = "Loading focalplanes from {}".format(fpdir) log.debug(msg) for root, dirs, files in os.walk(fpdir): for f in files: fpmat = fppat.match(f) if fpmat is not None: dt = if dt not in fpraw: fpraw[dt] = dict() fpraw[dt]["fp"] = os.path.join(root, f) continue stmat = stpat.match(f) if stmat is not None: dt = if dt not in fpraw: fpraw[dt] = dict() fpraw[dt]["st"] = os.path.join(root, f) continue exmat = expat.match(f) if exmat is not None: dt = if dt not in fpraw: fpraw[dt] = dict() fpraw[dt]["ex"] = os.path.join(root, f) break # Check that we have all 3 files needed for each timestamp for ts, files in fpraw.items(): for key in ["fp", "st", "ex"]: if key not in files: msg = "Focalplane state for time {} is missing one of \ the 3 required files (focalplane, state, exclusion)".format( ts ) raise RuntimeError(msg) # Now load the files for each time into our cached global variable. _focalplane = list() for ts in sorted(fpraw.keys()): dt = None file_dt = None try: file_dt = datetime.strptime(ts, "%Y-%m-%dT%H:%M:%S%z") dt = file_dt except ValueError: # This is an old file with implicit UTC times (no offset) # Load it as-is and then set the time zone. file_dt = datetime.strptime(ts, "%Y-%m-%dT%H:%M:%S") dt = file_dt.replace(tzinfo=timezone.utc) _focalplane.append(( dt, file_dt, { "fp_file": fpraw[ts]["fp"], "st_file": fpraw[ts]["st"], "ex_file": fpraw[ts]["ex"], "fp_data": None, "st_data": None, "ex_data": None, } ))
[docs]def get_focalplane_dates(): """Returns the dates of new focalplane definitions. There are two levels of time-dependent changes within the focalplane. First are the focalplane definitions, defined by a set of files ``$DESIMODEL/data/focalplane/{desi-exclusion,desi-focalplane,desi-state}_DATE.``. Those are the dates returned by this function, as a list of datetime objects. The second level is that, within the "state" table, there are changes to the states of individual positioners. (These are the dates returned when ``get_time_range=True`` is set in ``load_focalplane()``.) Returns ------- dates : :class:`list` of :class:`datetime` The dates when the focalplane changed. """ ensure_focalplane_loaded() global _focalplane return [dt for dt,fdt,v in _focalplane]
[docs]def load_focalplane(time=None, get_time_range=False): """Load the focalplane state that is valid for the given time. Parameters ---------- time : :class:`~datetime.datetime` The time to query with explicit timezone. Default to current time (:meth:``) with timezone UTC. Returns ------- :class:`tuple` A tuple of (FP layout, exclusion polygons, state, time string). The FP layout is a Table. The exclusion polygons are a dictionary indexed by names that are referenced in the state. The state is a Table. The time string is the resulting UTC ISO format time string for the creation date of the FP model. If get_time_range=True, returns two additional values: time_low and time_high, both datetime objects giving the range of dates over which this description of the focal plane is valid. `time_high` may be None, indicating that there is no later known hardware state. In particular, these dates refer to the `state` of the positioners, which are more fine-grained than the `fp` and `exclusion` objects. """ # Time range over which this FP model is valid. time_lo = time_hi = None if time is None: time = elif time.tzinfo is None: msg = "Requested focalplane time '{}' is not timezone-aware. Assuming UTC.".format(time) log.warning(msg) time = time.replace(tzinfo=timezone.utc) # Convert requested time to UTC time = time.astimezone(tz=timezone.utc) # Load the global _focalplane cache if required ensure_focalplane_loaded() # Search the list of states for the most recent time that is before our # requested time. There should not be too many different states, or else # we are using the wrong format for storing these. Therefore, a linear # search should be fast enough. focalplane_props = None file_tmstr = None for dt, file_dt, props in _focalplane: if time >= dt: focalplane_props = props try: file_tmstr = file_dt.isoformat(timespec="seconds") except TypeError: # This must be python < 3.6, with no timespec option. # Since the focalplane time is read from the file name without # microseconds, the microseconds should be zero and so the # default return string will be correct. file_tmstr = file_dt.isoformat() time_lo = dt else: time_hi = dt break if focalplane_props is None: msg = "Cannot find focalplane for time {}".format(time) raise RuntimeError(msg) # Load the data if this is the first time working with this focalplane if focalplane_props["fp_data"] is None: focalplane_props["fp_data"] = focalplane_props["fp_file"], format="ascii.ecsv" ) focalplane_props["st_data"] = focalplane_props["st_file"], format="ascii.ecsv" ) if (focalplane_props['ex_file'].endswith('.json') or focalplane_props['ex_file'].endswith('.json.gz')): import json loadroutine = json.load else: loadroutine = yaml.safe_load try: # First try to load uncompressed with open(focalplane_props["ex_file"], "r") as f: focalplane_props["ex_data"] = loadroutine(f) except: # Must be gzipped with["ex_file"], "rb") as f: focalplane_props["ex_data"] = loadroutine(f) # Now "replay" the state up to our requested time. st_data = focalplane_props["st_data"] locstate = dict() # Believe it or not, parsing strings into dates takes a significant fraction of the compute # time here, so cache them! parsed_dates = {} for row in range(len(st_data)): datestr = st_data[row]["TIME"] tm = parsed_dates.get(datestr) if tm is None: # not cached, parse it! try: tm = datetime.strptime(datestr, "%Y-%m-%dT%H:%M:%S%z") except ValueError: # Old format with implicit UTC timezone tm = datetime.strptime(datestr, "%Y-%m-%dT%H:%M:%S") tm = tm.replace(tzinfo=timezone.utc) parsed_dates[datestr] = tm if tm <= time: loc = st_data[row]["LOCATION"] locstate[loc] = st_data[row] time_lo = tm else: time_hi = tm break rows = list() for loc in sorted(locstate.keys()): rows.append(locstate[loc]) state_data = Table(rows=rows, names=st_data.colnames) state_data.remove_column("TIME") rtn = ( focalplane_props["fp_data"], focalplane_props["ex_data"], state_data, file_tmstr ) if get_time_range: return rtn + (time_lo, time_hi) return rtn
[docs]def reset_cache(): """Reset I/O cache.""" global _thru, _psf, _params, _gfa, _fiberpos, _tiles, _platescale, _focalplane _thru = dict() _psf = dict() _params = None _gfa = None _fiberpos = None _tiles = dict() _platescale = None _focalplane = None
[docs]def load_target_info(): """Loads data/targets/targets.yaml and returns the nested dictionary. This is primarily syntactic sugar to avoid end users constructing paths and filenames by hand (which *e.g.* broke when targets.dat was renamed to targets.yaml). Returns ------- :class:`dict` The dictionary read from the YAML file. """ targetsfile = findfile("targets/targets.yaml") if not os.path.exists(targetsfile): targetsfile = findfile("targets/targets.dat") with open(targetsfile) as fx: data = yaml.safe_load(fx) return data
[docs]def load_pixweight(nside, pixmap=None): """Loads ``$DESIMODEL/data/footprint/desi-healpix-weights.fits``. Parameters ---------- nside : :class:`int` After loading, the array will be resampled to the passed HEALPix `nside`. pixmap : :class:``, optional Input pixel weight map, already read from a weights file. Returns ------- Weight HEALPix weight map for the DESI footprint at the requested `nside`. """ import healpy as hp if pixmap is not None: log.debug("Using input pixel weight map of length {}.".format(len(pixmap))) else: # ADM read in the standard pixel weights file pixfile = findfile("footprint/desi-healpix-weights.fits") with as hdulist: pixmap = hdulist[0].data # ADM determine the file's nside, and flag a warning if the passed nside exceeds it npix = len(pixmap) truenside = hp.npix2nside(len(pixmap)) if truenside < nside: log.warning( "downsampling is fuzzy...Passed nside={}, but file {} is stored at nside={}".format( nside, pixfile, truenside ) ) # ADM resample the map return hp.pixelfunc.ud_grade(pixmap, nside, order_in="NESTED", order_out="NESTED")
[docs]def findfile(filename, surveyops=False): """Return full path to data file ``$DESIMODEL/data/filename``. Parameters ---------- filename : :class:`str` Name of the file, relative to the desimodel data directory. surveyops : :class:`bool` If ``True`` then find the relevant path for the $DESI_SURVEYOPS directory rather than the $DESIMODEL directory. Returns ------- :class:`str` The full path. Notes ----- This is a precursor for a potential future refactor where desimodel data would be installed with the package and :envvar:`DESIMODEL` would become an optional override. """ return os.path.join(datadir(surveyops), filename)
[docs]def datadir(surveyops=False): """Returns location to desimodel data. Parameters ---------- surveyops : :class:`bool` If ``True`` then find the relevant path for the $DESI_SURVEYOPS directory rather than the $DESIMODEL directory. Notes ----- If `surveyops`==``False`` and :envvar:`DESIMODEL` is set, then $DESIMODEL overrides data installed with the package. """ if surveyops: if "DESI_SURVEYOPS" in os.environ: surveyops = os.environ["DESI_SURVEYOPS"] # ADM test whether surveyops directory was checked out to trunk. if os.path.isdir(os.path.join(surveyops, "trunk", "ops")): surveyops = os.path.join(surveyops, "trunk") return os.path.abspath(os.path.join(surveyops, "ops")) # ADM raise a custom exception if $DESI_SURVEYOPS is not set. else: raise MissingEnvVar(f"$DESI_SURVEYOPS is not set") else: if "DESIMODEL" in os.environ: return os.path.abspath(os.path.join(os.environ["DESIMODEL"], "data")) else: import pkg_resources return pkg_resources.resource_filename("desimodel", "data")