"""
desisim.util
============
Utility functions for desisim. These may belong elsewhere?
"""
import numpy as np
[docs]
def hadec2airmass(ha, dec, latitude=31.96397222):
"""
Return airmass of an observation at hour angle `ha` and declination `dec`
Args:
ha (float): hour angle in degrees
dec (float): declination in degrees
Options:
latitude (float): telescope latitude; defaults to Kitt Peak
Returns: airmass (float)
Code adapted from desisurvey.util
"""
# convert inputs from degrees to radians
ha = np.radians(ha)
dec = np.radians(dec)
latitude = np.radians(latitude)
# cos(zenith)
cosZ = (np.sin(dec) * np.sin(latitude) +
np.cos(dec) * np.cos(latitude) * np.cos(ha))
# cos(zenith) -> airmass, following Rozenberg 1966
# see https://en.wikipedia.org/wiki/Air_mass_(astronomy)#Interpolative_formulas
airmass = 1. / (cosZ + 0.025 * np.exp(-11 * cosZ))
# ensure that rounding doesn't result in airmass<1
airmass = np.clip(airmass, 1.0, None)
return airmass
[docs]
def spline_medfilt2d(image, kernel_size=201):
'''
Returns a 2D spline interpolation of a median filtered input image
'''
if 3*kernel_size > min(image.shape):
raise ValueError(
'kernel_size {} must be < min image shape {}//3'.format(
kernel_size, min(image.shape)))
from scipy.interpolate import RectBivariateSpline
n = kernel_size // 2
xx = np.arange(n, image.shape[1], kernel_size)
yy = np.arange(n, image.shape[0], kernel_size)
zz = np.zeros((len(yy), len(xx)))
for i,x in enumerate(xx):
for j,y in enumerate(yy):
xy = np.s_[y-n:y+n+1, x-n:x+n+1]
zz[j,i] = np.median(image[xy])
#- adjust spline order for small test data
kx = min(3, len(xx)-1)
ky = min(3, len(yy)-1)
s = RectBivariateSpline(xx, yy, zz, kx=kx, ky=ky)
background = s(np.arange(image.shape[0]), np.arange(image.shape[1]))
return background
[docs]
def medxbin(x,y,binsize,minpts=20,xmin=None,xmax=None):
"""
Compute the median (and other statistics) in fixed bins along the x-axis.
"""
import numpy as np
from scipy import ptp
# Need an exception if there are fewer than three arguments.
if xmin==None:
xmin = x.min()
if xmax==None:
xmax = x.max()
#print(xmin,xmax)
nbin = int(ptp(x)/binsize)
bins = np.linspace(xmin,xmax,nbin)
idx = np.digitize(x,bins)
#print(nbin, bins, xmin, xmax)
stats = np.zeros(nbin,[('median','f8'),('sigma','f8'),('iqr','f8')])
for kk in np.arange(nbin):
npts = len(y[idx==kk])
if npts>minpts:
stats['median'][kk] = np.median(y[idx==kk])
stats['sigma'][kk] = np.std(y[idx==kk])
stats['iqr'][kk] = np.subtract(*np.percentile(y[idx==kk],[75, 25]))
# Remove bins with too few points.
good = np.nonzero(stats['median'])
stats = stats[good]
return bins[good], stats
#- TODO: move to desiutil
[docs]
def dateobs2night(dateobs):
'''
Convert UTC dateobs to KPNO YEARMMDD night string
Args:
dateobs:
float -> interpret as MJD
str -> interpret as ISO 8601 YEAR-MM-DDThh:mm:ss.s string
astropy.time.Time -> UTC
python datetime.datetime -> UTC
TODO: consider adding format option to pass to astropy.time.Time without
otherwise questioning the dateobs format
'''
# See pixsim.py
from desiutil.iers import freeze_iers
freeze_iers()
import astropy.time
import datetime
if isinstance(dateobs, float):
dateobs = astropy.time.Time(dateobs, format='mjd')
elif isinstance(dateobs, datetime.datetime):
dateobs = astropy.time.Time(dateobs, format='datetime')
elif isinstance(dateobs, str):
dateobs = astropy.time.Time(dateobs, format='isot')
elif not isinstance(dateobs, astropy.time.Time):
raise ValueError('dateobs must be float, str, datetime, or astropy time object')
import astropy.units as u
kpno_time = dateobs - 7*u.hour
#- "night" rolls over at local noon, not midnight, so subtract another 12 hours
yearmmdd = (kpno_time - 12*u.hour).isot[0:10].replace('-', '')
assert len(yearmmdd) == 8
return yearmmdd