Spherical harmonic maps¶
- class starry.Map
The default
starry
map class.This class handles light curves and phase curves of objects in emitted light.
Note
Instantiate this class by calling
starry.Map()
withydeg > 0
and bothrv
andreflected
set to False.- property N
Total number of map coefficients. Read-only
This is equal to \(N_\mathrm{y} + N_\mathrm{u} + N_\mathrm{f}\).
- property Nf
Number of spherical harmonic coefficients in the filter. Read-only
This is equal to \((f_\mathrm{deg} + 1)^2\).
- property Nu
Number of limb darkening coefficients, including \(u_0\). Read-only
This is equal to \(u_\mathrm{deg} + 1\).
- property Ny
Number of spherical harmonic coefficients. Read-only
This is equal to \((y_\mathrm{deg} + 1)^2\).
- add_spot(amp=None, intensity=None, relative=True, sigma=0.1, lat=0.0, lon=0.0)
Add the expansion of a gaussian spot to the map.
This function adds a spot whose functional form is the spherical harmonic expansion of a gaussian in the quantity \(\cos\Delta\theta\), where \(\Delta\theta\) is the angular separation between the center of the spot and another point on the surface. The spot brightness is controlled by either the parameter
amp
, defined as the fractional change in the total luminosity of the object due to the spot, or the parameterintensity
, defined as the fractional change in the intensity at the center of the spot.- Parameters
amp (scalar or vector, optional) – The amplitude of the spot. This is equal to the fractional change in the luminosity of the map due to the spot. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Default is None. Either
amp
orintensity
must be given.intensity (scalar or vector, optional) – The intensity of the spot. This is equal to the fractional change in the intensity of the map at the center of the spot. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Default is None. Either
amp
orintensity
must be given.relative (bool, optional) – If True, computes the spot expansion assuming the fractional amp or intensity change is relative to the current map amplitude/intensity. If False, computes the spot expansion assuming the fractional change is relative to the original map amplitude/intensity (i.e., that of a featureless map). Defaults to True. Note that if True, adding two spots with the same values of amp or intensity will generally result in different intensities at their centers, since the first spot will have changed the map intensity everywhere! Defaults to True.
sigma (scalar, optional) – The standard deviation of the gaussian. Defaults to 0.1.
lat (scalar, optional) – The latitude of the spot in units of
angle_unit
. Defaults to 0.0.lon (scalar, optional) – The longitude of the spot in units of
angle_unit
. Defaults to 0.0.
Warning
Method
add_spot
is deprecated as of version1.1
and will be removed in the future. Please use methodspot
instead.
- property amp
The overall amplitude of the map in arbitrary units. This factor multiplies the intensity and the flux and is thus proportional to the luminosity of the object. For multi-wavelength maps, this is a vector corresponding to the amplitude of each wavelength bin. For reflected light maps, this is the average spherical albedo of the body.
- property angle_unit
An
astropy.units
quantity defining the angle unit for this map.
- property deg
Total degree of the map. Read-only
This is equal to \(y_\mathrm{deg} + u_\mathrm{deg} + f_\mathrm{deg}\).
- design_matrix(**kwargs)
Compute and return the light curve design matrix \(A\).
This matrix is used to compute the flux \(f\) from a vector of spherical harmonic coefficients \(y\) and the map amplitude \(a\): \(f = a A y\).
- Parameters
xo (scalar or vector, optional) – x coordinate of the occultor relative to this body in units of this body’s radius.
yo (scalar or vector, optional) – y coordinate of the occultor relative to this body in units of this body’s radius.
zo (scalar or vector, optional) – z coordinate of the occultor relative to this body in units of this body’s radius.
ro (scalar, optional) – Radius of the occultor in units of this body’s radius.
theta (scalar or vector, optional) – Angular phase of the body in units of
angle_unit
.
- draw()
Draw a map from the posterior distribution.
This method draws a random map from the posterior distribution and sets the
y
map vector andamp
map amplitude accordingly. Users should callsolve()
to enable this attribute.
- property fdeg
Degree of the multiplicative filter. Read-only
- flux(**kwargs)
Compute and return the light curve.
- Parameters
xo (scalar or vector, optional) – x coordinate of the occultor relative to this body in units of this body’s radius.
yo (scalar or vector, optional) – y coordinate of the occultor relative to this body in units of this body’s radius.
zo (scalar or vector, optional) – z coordinate of the occultor relative to this body in units of this body’s radius.
ro (scalar, optional) – Radius of the occultor in units of this body’s radius.
theta (scalar or vector, optional) – Angular phase of the body in units of
angle_unit
.integrated (bool, optional) – If True, dots the flux with the amplitude. Default False, in which case this returns a 2d array (wavelength-dependent maps only).
- get_latlon_grid(res=300, projection='ortho')
Return the latitude/longitude grid corresponding to the result of a call to
render()
.- Parameters
res (int, optional) – The resolution of the map in pixels on a side. Defaults to 300.
projection (string, optional) – The map projection. Accepted values are
ortho
, corresponding to an orthographic projection (as seen on the sky),rect
, corresponding to an equirectangular latitude-longitude projection, andmoll
, corresponding to a Mollweide equal-area projection. Defaults toortho
.
- get_pixel_transforms(oversample=2, lam=1e-06, eps=1e-06)
Return several linear operators for pixel transformations.
- Parameters
oversample (int) – Factor by which to oversample the pixelization grid. Default 2.
lam (float) – Regularization parameter for the inverse pixel transform. Default 1e-6.
eps (float) – Regularization parameter for the derivative transforms. Default 1e-6.
- Returns
The tuple (lat, lon, Y2P, P2Y, Dx, Dy).
The transforms returned by this method can be used to easily convert back and forth between spherical harmonic coefficients and intensities on a discrete pixelized grid. Projections onto pixels are performed on an equal-area Mollweide grid, so these transforms are useful for applying priors on the pixel intensities, for instance.
The lat and lon arrays correspond to the latitude and longitude of each of the points used in the transform (in units of angle_unit).
The Y2P matrix is an operator that transforms from spherical harmonic coefficients y to pixels p on a Mollweide grid:
p = Y2P @ y
The P2Y matrix is the (pseudo-)inverse of that operator:
y = P2Y @ p
Finally, the Dx and Dy operators transform a pixel representation of the map p to the derivative of p with respect to longitude and latitude, respectively:
dpdlon = Dx @ p dpdlat = Dy @ p
By combining these operators, one can differentiate the spherical harmonic expansion with respect to latitude and longitude, if desired:
dydlon = P2Y @ Dx @ Y2P @ y dydlat = P2Y @ Dy @ Y2P @ y
These derivatives could be useful for implementing total-variation-reducing regularization, for instance.
Warning
This is an experimental feature.
- property inc
The inclination of the rotation axis in units of
angle_unit
.
- intensity(lat=0, lon=0, **kwargs)
Compute and return the intensity of the map.
- Parameters
lat (scalar or vector, optional) – latitude at which to evaluate the intensity in units of
angle_unit
.lon (scalar or vector, optional) – longitude at which to evaluate the intensity in units of
angle_unit`
.theta (scalar, optional) – For differentially rotating maps only, the angular phase at which to evaluate the intensity. Default 0.
limbdarken (bool, optional) – Apply limb darkening (only if
udeg
> 0)? Default True.
- intensity_design_matrix(lat=0, lon=0)
Compute and return the pixelization matrix
P
.This matrix is used to compute the intensity \(I\) from a vector of spherical harmonic coefficients \(y\) and the map amplitude \(a\): \(I = a P y\).
- Parameters
lat (scalar or vector, optional) – latitude at which to evaluate the design matrix in units of
angle_unit
.lon (scalar or vector, optional) – longitude at which to evaluate the design matrix in units of
angle_unit
.
Note
This method ignores any filters (such as limb darkening or velocity weighting) and illumination (for reflected light maps).
- limbdark_is_physical()
Check whether the limb darkening profile (if any) is physical.
This method uses Sturm’s theorem to ensure that the limb darkening intensity is positive everywhere and decreases monotonically toward the limb.
- Returns
Whether or not the limb darkening profile is physical.
- Return type
bool
- lnlike(*, design_matrix=None, woodbury=True, **kwargs)
Returns the log marginal likelihood of the data given a design matrix.
This method computes the marginal likelihood (marginalized over the spherical harmonic coefficients) given a light curve and its covariance (set via the
set_data()
method) and a Gaussian prior on the spherical harmonic coefficients (set via theset_prior()
method).- Parameters
design_matrix (matrix, optional) – The flux design matrix, the quantity returned by
design_matrix()
. Default is None, in which case this is computed based onkwargs
.woodbury (bool, optional) – Solve the linear problem using the Woodbury identity? Default is True. The Woodbury identity is used to speed up matrix operations in the case that the number of data points is much larger than the number of spherical harmonic coefficients. In this limit, it can speed up the code by more than an order of magnitude. Keep in mind that the numerical stability of the Woodbury identity is not great, so if you’re getting strange results try disabling this. It’s also a good idea to disable this in the limit of few data points and large spherical harmonic degree.
kwargs (optional) – Keyword arguments to be passed directly to
design_matrix()
, if a design matrix is not provided.
- Returns
The log marginal likelihood, a scalar.
- load(image, extent=(- 180, 180, - 90, 90), smoothing=None, fac=1.0, eps=1e-12, force_psd=False, **kwargs)
Load an image or ndarray.
This routine performs a simple spherical harmonic transform (SHT) to compute the spherical harmonic expansion corresponding to an input image file or
numpy
array on a lat-lon grid. The resulting coefficients are ingested into the map.- Parameters
image – A path to an image PNG file or a two-dimensional
numpy
array on a latitude-longitude grid.extent (tuple, optional) – The lat-lon values corresponding to the edges of the image in degrees,
(lat0, lat1, lon0, lon1)
. Default is(-180, 180, -90, 90)
.smoothing (float, optional) – Gaussian smoothing strength. Increase this value to suppress ringing or explicitly set to zero to disable smoothing. Default is
1/self.ydeg
.fac (float, optional) – Factor by which to oversample the image when applying the SHT. Default is
1.0
. Increase this number for higher fidelity (at the expense of increased computational time).eps (float, optional) – Regularization strength for the spherical harmonic transform. Default is
1e-12
.force_psd (bool, optional) – Force the map to be positive semi-definite? Default is False.
kwargs (optional) – Any other kwargs passed directly to
minimize()
(only ifpsd
is True).
- minimize(oversample=1, ntries=1, bounds=None, return_info=False)
Find the global (optionally local) minimum of the map intensity.
- Parameters
oversample (int) – Factor by which to oversample the initial grid on which the brute force search is performed. Default 1.
ntries (int) – Number of times the nonlinear minimizer is called. Default 1.
return_info (bool) – Return the info from the minimization call? Default is False.
bounds (tuple) – Return map minimum in a certain latitude/longitude range, for example bounds=((0, 90), (0, 180)). Default None.
- Returns
A tuple of the latitude, longitude, and the value of the intensity at the minimum. If
return_info
is True, also returns the detailed solver information.
- property nw
Number of wavelength bins. Read-only
- property obl
The obliquity of the rotation axis in units of
angle_unit
.
- remove_prior()
Remove the prior on the map coefficients.
- render(res=300, projection='ortho', theta=0.0)
Compute and return the intensity of the map on a grid.
Returns an image of shape
(res, res)
, unlesstheta
is a vector, in which case returns an array of shape(nframes, res, res)
, wherenframes
is the number of values oftheta
. However, if this is a spectral map,nframes
is the number of wavelength bins andtheta
must be a scalar.Note
Users can obtain the latitudes and longitudes corresponding to each point in the rendered image by calling
get_latlon_grid()
.- Parameters
res (int, optional) – The resolution of the map in pixels on a side. Defaults to 300.
projection (string, optional) – The map projection. Accepted values are
ortho
, corresponding to an orthographic projection (as seen on the sky),rect
, corresponding to an equirectangular latitude-longitude projection, andmoll
, corresponding to a Mollweide equal-area projection. Defaults toortho
.theta (scalar or vector, optional) – The map rotation phase in units of
angle_unit
. If this is a vector, an animation is generated. Defaults to0.0
.
- reset(**kwargs)
Reset all map coefficients and attributes.
Note
Does not reset custom unit settings.
- rotate(axis, theta)
Rotate the current map vector an angle
theta
aboutaxis
.- Parameters
axis (vector) – The axis about which to rotate the map.
theta (scalar) – The angle of (counter-clockwise) rotation.
- set_data(flux, C=None, cho_C=None)
Set the data vector and covariance matrix.
This method is required by the
solve()
method, which analytically computes the posterior over surface maps given a dataset and a prior, provided both are described as multivariate Gaussians.- Parameters
flux (vector) – The observed light curve.
C (scalar, vector, or matrix) – The data covariance. This may be a scalar, in which case the noise is assumed to be homoscedastic, a vector, in which case the covariance is assumed to be diagonal, or a matrix specifying the full covariance of the dataset. Default is None. Either C or cho_C must be provided.
cho_C (matrix) – The lower Cholesky factorization of the data covariance matrix. Defaults to None. Either C or cho_C must be provided.
- set_prior(*, mu=None, L=None, cho_L=None)
Set the prior mean and covariance of the spherical harmonic coefficients.
This method is required by the
solve()
method, which analytically computes the posterior over surface maps given a dataset and a prior, provided both are described as multivariate Gaussians.Note that the prior is placed on the amplitude-weighted coefficients, i.e., the quantity
x = map.amp * map.y
. Because the first spherical harmonic coefficient is fixed at unity,x[0]
is the amplitude of the map. The actual spherical harmonic coefficients are given byx / map.amp
.This convention allows one to linearly fit for an arbitrary map normalization at the same time as the spherical harmonic coefficients, while ensuring the
starry
requirement that the coefficient of the \(Y_{0,0}\) harmonic is always unity.- Parameters
mu (scalar or vector) – The prior mean on the amplitude-weighted spherical harmonic coefficients. Default is 1.0 for the first term and zero for the remaining terms. If this is a vector, it must have length equal to
Ny
.L (scalar, vector, or matrix) – The prior covariance. This may be a scalar, in which case the covariance is assumed to be homoscedastic, a vector, in which case the covariance is assumed to be diagonal, or a matrix specifying the full prior covariance. Default is None. Either L or cho_L must be provided.
cho_L (matrix) – The lower Cholesky factorization of the prior covariance matrix. Defaults to None. Either L or cho_L must be provided.
- show(**kwargs)
Display an image of the map, with optional animation. See the docstring of
render()
for more details and additional keywords accepted by this method.- Parameters
ax (optional) – A matplotlib axis instance to use. Default is to create a new figure.
cmap (string or colormap instance, optional) – The matplotlib colormap to use. Defaults to
plasma
.figsize (tuple, optional) – Figure size in inches. Default is (3, 3) for orthographic maps and (7, 3.5) for rectangular maps.
projection (string, optional) – The map projection. Accepted values are
ortho
, corresponding to an orthographic projection (as seen on the sky),rect
, corresponding to an equirectangular latitude-longitude projection, andmoll
, corresponding to a Mollweide equal-area projection. Defaults toortho
.colorbar (bool, optional) – Display a colorbar? Default is False.
grid (bool, optional) – Show latitude/longitude grid lines? Defaults to True.
interval (int, optional) – Interval between frames in milliseconds (animated maps only). Defaults to 75.
file (string, optional) – The file name (including the extension) to save the figure or animation to. Defaults to None.
html5_video (bool, optional) – If rendering in a Jupyter notebook, display as an HTML5 video? Default is True. If False, displays the animation using Javascript (file size will be larger.)
dpi (int, optional) – Image resolution in dots per square inch. Defaults to the value defined in
matplotlib.rcParams
.bitrate (int, optional) – Bitrate in kbps (animations only). Defaults to the value defined in
matplotlib.rcParams
.norm (optional) – The color normalization passed to
matplotlib.pyplot.imshow
, an instance ofmatplotlib.colors.Normalize
. Can be used to pass in minimum and maximum values. Default is None.illuminate (bool, optional) – Illuminate the map (reflected light maps only)? Default True. If False, shows the albedo surface map.
screen (bool, optional) – Apply the illumination as a black-and-white alpha screen (reflected light maps only)? Default True. If False, a single colormap is used to plot the visible intensity.
Note
Pure limb-darkened maps do not accept a
projection
keyword.Note
If calling this method on an instance of
Map
created within apymc3.Model()
, you may specify apoint
keyword with the model point at which to evaluate the map. This method also accepts amodel
keyword, although this is inferred automatically if called from within apymc3.Model()
context. If no point is provided, attempts to evaluate the map atmodel.test_point
and raises a warning.
- sht_matrix(inverse=False, return_grid=False, smoothing=None, oversample=2, lam=1e-06)
Return the Spherical Harmonic Transform (SHT) matrix.
This matrix dots into a vector of pixel intensities defined on a set of latitude-longitude points, resulting in a vector of spherical harmonic coefficients that best approximates the image.
This can be useful for transforming between the spherical harmonic and pixel representations of the surface map, such as when specifying priors during inference or optimization. For example, one application is when one wishes to sample in the pixel intensities
p
instead of the spherical harmonic coefficientsy
. In this case, one must compute the spherical harmonic coefficientsy
from the vector of pixel intensitiesp
so thatstarry
can compute the model for the flux:Note
This method is a thin wrapper of the
get_pixel_transforms
method.- Parameters
inverse (bool, optional) – which transforms from spherical harmonic coefficients to pixels. Default is False.
return_grid (bool, optional) – shape
(npix, 2)
corresponding to the latitude-longitude points (in units ofangle_unit
) on which the SHT is evaluated. Default is False.smoothing (float, optional) – Gaussian smoothing strength. Increase this value to suppress ringing (forward SHT only) or explicitly set to zero to disable smoothing. Default is
2/self.ydeg
.oversample (int, optional) – Factor by which to oversample the pixelization grid. Default 2.
lam (float, optional) – Regularization parameter for the inverse pixel transform. Default 1e-6.
- Returns
A matrix of shape (
Ny
,npix
) or (npix
,Ny
) (ifinverse
is True), wherenpix
is the number of pixels on the grid. This number is determined from the degree of the map and theoversample
keyword. Ifreturn_grid
is True, also returns the latitude-longitude points (in units ofangle_unit
) on which the SHT is evaluated, a matrix of shape(npix, 2)
.
- property solution
The posterior probability distribution for the map.
This is a tuple containing the mean and lower Cholesky factorization of the covariance of the amplitude-weighted spherical harmonic coefficient vector, obtained by solving the regularized least-squares problem via the
solve()
method.Note that to obtain the actual covariance matrix from the lower Cholesky factorization \(L\), simply compute \(L L^\top\).
Note also that this is the posterior for the amplitude-weighted map vector. Under this convention, the map amplitude is equal to the first term of the vector and the spherical harmonic coefficients are equal to the vector normalized by the first term.
- solve(*, design_matrix=None, **kwargs)
Solve the linear least-squares problem for the posterior over maps.
This method solves the generalized least squares problem given a light curve and its covariance (set via the
set_data()
method) and a Gaussian prior on the spherical harmonic coefficients (set via theset_prior()
method). The map amplitude and coefficients are set to the maximum a posteriori (MAP) solution.- Parameters
design_matrix (matrix, optional) – The flux design matrix, the quantity returned by
design_matrix()
. Default is None, in which case this is computed based onkwargs
.kwargs (optional) – Keyword arguments to be passed directly to
design_matrix()
, if a design matrix is not provided.
- Returns
A tuple containing the posterior mean for the amplitude-weighted spherical harmonic coefficients (a vector) and the Cholesky factorization of the posterior covariance (a lower triangular matrix).
Note
Users may call
draw()
to draw from the posterior after calling this method.
- spot(*, contrast=1.0, radius=None, lat=0.0, lon=0.0, **kwargs)
Add the expansion of a circular spot to the map.
This function adds a spot whose functional form is a top hat in \(\Delta\theta\), the angular separation between the center of the spot and another point on the surface. The spot intensity is controlled by the parameter
contrast
, defined as the fractional change in the intensity at the center of the spot.- Parameters
contrast (scalar or vector, optional) – The contrast of the spot. This is equal to the fractional change in the intensity of the map at the center of the spot relative to the baseline intensity of an unspotted map. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Positive values of the contrast result in dark spots; negative values result in bright spots. Default is
1.0
, corresponding to a spot with central intensity close to zero.radius (scalar, optional) – The angular radius of the spot in units of
angle_unit
. Defaults to20.0
degrees.lat (scalar, optional) – The latitude of the spot in units of
angle_unit
. Defaults to0.0
.lon (scalar, optional) – The longitude of the spot in units of
angle_unit
. Defaults to0.0
.
Note
Keep in mind that things are normalized in
starry
such that the disk-integrated flux (not the intensity!) of an unspotted body is unity. The default intensity of an unspotted map is1.0 / np.pi
everywhere (this ensures the integral over the unit disk is unity). So when you instantiate a map and add a spot of contrastc
, you’ll see that the intensity at the center is actually(1 - c) / np.pi
. This is expected behavior, since that’s a factor of1 - c
smaller than the baseline intensity.Note
This function computes the spherical harmonic expansion of a circular spot with uniform contrast. At finite spherical harmonic degree, this will return an approximation that may be subject to ringing. Users can control the amount of ringing and the smoothness of the spot profile (see below). In general, however, at a given spherical harmonic degree
ydeg
, there is always minimum spot radius that can be modeled well. Forydeg = 15
, for instance, that radius is about10
degrees. Attempting to add a spot smaller than this will in general result in a large amount of ringing and a smaller contrast than desired.There are a few additional under-the-hood keywords that control the behavior of the spot expansion. These are
- Parameters
spot_pts (int, optional) – The number of points in the expansion of the (1-dimensional) spot profile. Default is
1000
.spot_eps (float, optional) – Regularization parameter in the expansion. Default is
1e-9
.spot_smoothing (float, optional) – Standard deviation of the Gaussian smoothing applied to the spot to suppress ringing (unitless). Default is
2.0 / self.ydeg
.spot_fac (float, optional) – Parameter controlling the smoothness of the spot profile. Increasing this parameter increases the steepness of the profile (which approaches a top hat as
spot_fac -> inf
). Decreasing it results in a smoother sigmoidal function. Default is300
. Changing this parameter is not recommended; changespot_smoothing
instead.
Note
These last four parameters are cached. That means that changing their value in a call to
spot
will result in all future calls tospot
“remembering” those settings, unless you change them back!
- property u
The vector of limb darkening coefficients. Read-only
To set this vector, index the map directly using one index:
map[n] = ...
wheren
is the degree of the limb darkening coefficient. This may be an integer or an array of integers. Slice notation may also be used.
- property udeg
Limb darkening degree. Read-only
- property wav
The wavelength(s) at which the flux is measured in units of
wav_unit
.
- property wav_unit
An
astropy.units
quantity defining the angle unit for this map.
- property y
The spherical harmonic coefficient vector. Read-only
To set this vector, index the map directly using two indices:
map[l, m] = ...
wherel
is the spherical harmonic degree andm
is the spherical harmonic order. These may be integers or arrays of integers. Slice notation may also be used.
- property ydeg
Spherical harmonic degree of the map. Read-only