PSF¶

class snappl.psf.PSF(x=None, y=None, band=None, observation_id=None, sca=None, _called_from_get_psf_object=False, image=None, seed=None, **kwargs)[source]¶

Bases: object

Wraps a PSF. All roman snpit photometry code will ideally only use PSF methods defined in this base class.

This is an abstract base class; it can do almost nothing itself. In practice, you need to instantiate a subclass. Do that by calling the class method PSF.get_psf_object.

Don’t call this or the constructor of a subclass directly, call PSF.get_psf_object().

See get_psf_object for parameter documentation.

_called_from_get_psf_object is used internally and should not be used outside this module, unless you know what you’re doing and intentionally mean to subvert the system.

Note:¶

If both are set, then band overrides image.band

Attributes Summary

`clip_size`
`stamp_size`	The size of the one side of a PSF image stamp at image resolution.
`x`
`y`

Methods Summary

`getImagePSF`([imagesampled])	Return a photutils.psf.ImagePSF model.
`get_psf_object`(psfclass[, x, y, band, ...])	Return a PSF object whose type is specified by psfclass.
`get_stamp`([x, y, x0, y0, flux])	Return a 2d numpy image of the PSF at the image resolution.

Attributes Documentation

clip_size¶

stamp_size¶: The size of the one side of a PSF image stamp at image resolution. Is always odd.

x¶

y¶

Methods Documentation

getImagePSF(imagesampled=True)[source]¶

Return a photutils.psf.ImagePSF model.

This is useful if you want to do, e.g., PSF photometry with photutils.

WARNING: at least with the photutilsImagePSF class, if you constructed it with peakx and peaky not at their defaults, then what you get back may not work as you expect. (It’s most usual to just leave peakx and peaky at their defaults; if you’re not doing that, then be very careful calling getImagePSF.)

Parameters:: imagesampled (bool, default) – By default, this getImagePSF a PSF model at image resolution. Set imagesampled to False, and it might return an oversampled PSF model; this will be class- and object-dependent.
Return type:: photutils.psf.ImagePSF

classmethod get_psf_object(psfclass, x=None, y=None, band=None, observation_id=None, sca=None, image=None, seed=None, **kwargs)[source]¶

Return a PSF object whose type is specified by psfclass.

Parameters:

psfclass (str) – The name of the class of PSF you want. Current options are: * photutilsImagePSF – a wrapper (sort of) around photutils.psf.ImagePSF * OversampledImagePSF – a PSF defined by an image which may be oversampled relative to the host image * YamlSerialized_OversampledImagePSF – OversampledImagePSF with a defined save formatr * A25ePSF – YamlSearialized_OversampledImagePSF from Aldoroty et al. 2025 * ou24PSF_slow – a PSF from galsim for OpenUniverse 2024 * ou24PSF – a PSF from galsim for OpenUniverse 2024 * STPSF – a PSF from STScI STPSF
x (float) –
The position on the host image that this is the PSF for. Usually you want these to have no fractional part (so x==floor(x) and y==floor(y)), meaning that you’ve evaluated the PSDF at the center of a pixel. (If you have PFSs that vary significantly on less than a pixel scale, you have such big problems that you probably shouldn’t even be trying to do astronomy.) Sometimes, but rarely. there is a use case for these values to have a on-zero fractional part.

Will default to something sane.

The exact definition of how this is used currently (sadly) depends a bit on the subclass. See the subclass constructor documentation for more details.
y (float) –
The position on the host image that this is the PSF for. Usually you want these to have no fractional part (so x==floor(x) and y==floor(y)), meaning that you’ve evaluated the PSDF at the center of a pixel. (If you have PFSs that vary significantly on less than a pixel scale, you have such big problems that you probably shouldn’t even be trying to do astronomy.) Sometimes, but rarely. there is a use case for these values to have a on-zero fractional part.

Will default to something sane.

The exact definition of how this is used currently (sadly) depends a bit on the subclass. See the subclass constructor documentation for more details.
band (str) – The Roman band this is a PSF for. (I.e. the band of the host image.) Ignored by some subclasses.
observation_id (str, default None) – Roman observation id. Ignored by many subclasses. (OU2024 uses this for its “pointing”.)
sca (int, default None) – Probably only relevant for the ou24PSF classes.
image (snappl.image.Image) –
Optional, The image that the psf is for. This is not used all subclasses, but is needed for some, currently the ou24 PSFs.

seed: int, default None
A random seed to pass to galsim.BaseDeviate for photonOps.
**kwargs (further keyword arguments passed to object constructor) – Specific subclasses may require or accept additional arguments. They will be documented in the subclasses’s constructor.
psfclass – The name of the class of the PSF to instantiate.
**kwargs

get_stamp(x=None, y=None, x0=None, y0=None, flux=1.0)[source]¶

Return a 2d numpy image of the PSF at the image resolution.

There are a distressing number of subtleties here, warranting an extended discussion.

INDEXING IMAGES¶

For discussion of pixel positions below, remember the conventions for astronomical arrays. Consider four things:

First thing to consider: in python, numpy arrays are 0-indexed. That is, if you have a 3-element numpy array named arr, the first element of the array is arr[0], the second arr[1], and the last arr[2]. Some other languages (e.g. FORTRAN) assume 1-indexed arrays. That is, the first element of FORTRAN array A is A[1], not A[0]. This matters for us because we are using some astronomical formats that have been around since everybody spoke Latin and everybody programmed in FORTRAN, so there are some legacy conventions left over. Some libraries (e.g. galsim) at least sometimes require you to specify array indexes (such as pixel positions) assuming 1-indexed arrays. Be very careful and read lots of documentation! If we’ve done it right, everything in snappl uses standard python numpy 0-based array indexes, so you will hopefully not become confused. What’s more more, the astropy.wcs.WCS class also uses the convention of 0-based arrays. (However, be careful, because astropy.wcs has an alternate interface that uses the other convention.)

Another place you will find 1-indexed arrays are in the WCSes defined in FITS headers, and in at least some FITS image display programs. If you use ds9 (the standard FITS image display program), and hover your pointer over the center of the lower-left pixel, you will notice that it tells you it’s at (1.0,1.0). This means that if you’re reading positions off of ds9, you always have to be careful to mentally convert when comparing to positions in your code! Likewise, if you try to manually apply the WCS transformation from a FITS header (doing the matrix multiplication yourself, rather than relying on a snappl or astropy library), you have to make sure you’re using 1-offset pixel coordinates. Generally, you will not have to worry about this; the WCS classes in snappl (just as in astropy) will internally take care of all these off-by-1 errors. As stated above, all snappl classes assume 0-based array indexing.

If you find yourself manually correcting for 1-offset pixel positions in your code, there’s a good chance you’re doing it wrong. snappl is supposed to take care of all of that.

Second thing to consider: following how numpy arrays are indexed, the lower-left pixel of an astronomical image is at x=0, y=0. Furthermore, by convention, the center of the lower-left pixel is at x=0.0, y=0.0. That means that for a 512×512 image, the whole array spans (-0.5,-0.5) to (511.5,511.5); the lower-left corner of the array, which is the lower-left corner of the lower-left pixel, is at (-0.5,-0.5).

Third thing to consider: because numpy arrays are (by default) stored in “row major” format, their indexing is backwards from what we might expect. That is, to get to pixel (x,y) on a numpy array image, you’d do image[y,x].

Fourth thing to consider: it follows that a pixel position whose fractional part is exactly 0.5 is right on the edge between two pixels. For example, the position x=0.5, y=0.5 is the corner between the four lower-left-most pixels on the image. If you want to ask for “the closest pixel center” in this case, there is an ambiguity, so we have to pick a convention; that convention is described below.

PSF CENTERING FOR get_stamp¶

If (x0, y0) are not given, the PSF will be centered as best possible on the stamp*†. So, if x ends in 0.8, it will be left of center, and if x ends in 0.2, it will be right of center. If the fractional part of x or y is exactly 0.5, there’s an ambituity as to where on the image you should place the stamp of the PSF. The position of the PSF on the returned stamp will always round down in this case. (The pixel on the image that corresponds to the center pixel on the stamp is at floor(x+0.5),floor(y+0.5), not round(x+0.5),round(y+0.5). Those two things are different, and round is not consistent. round(i.5) will round up if i is odd, but down if i is even. This makes it very difficult to understand where your PSF is; by using floor(x+0.5), we get consistent results regardless of whether the integer part of x is even or odd.)

For further discusison of centering, see the discusison of the (x0, y0) parameters below.

“The PSF will be centered as best possible on the stamp”: this is only true if the PSF itself is intrinsically centered. It’s possible that some subclasses will have non-intrinsically-centered PSFs. See the documentation on the __init__ and get_stamp methods of those subclasses (e.g. OversampledImagePSF and photutilsImagePSF) to make sure you understand how each subclass handles those cases. In all cases, get_stamp should return stamps that are consistent with the description in this docstring. If a subclass does something different, that subclass is broken.

† “Centered” is obvious when a PSF is perfectly radially: symmetric: the center of the PSF is its peak, or mode. If the PSF is not radially symmetric, then this becomes potentially ambiguous. The “center” of the PSF really becomes a “fiducial point”, and cannot be assumed to be the centroid or mode of the PSF (and the centroid and mode may well be different in this case). Hopefully it’s somewhere close. If you use consistent PSFs, then relative positions should be realiable. That is, if you do a PSF fit to an image to find positions of stars, and use the PSF positions of those stars with a WCS to find ra and dec, this will only work if you used the same PSFs to find the standard stars you used to solve for the WCS! For most of this discussion, for simplicitly, we’ll be assuming a radially symmetric PSF so that “peak” and “center” and “fiducial point” all mean the same thing.

param x:

Position on the image of the center of the psf. If not given, defaults to something sensible that was defined when the object was constructed. If you want to do sub-pixel shifts, then the fractional part of x will (usually) not be 0.

type x:

floats

param y:

Position on the image of the center of the psf. If not given, defaults to something sensible that was defined when the object was constructed. If you want to do sub-pixel shifts, then the fractional part of x will (usually) not be 0.

type y:

floats

param x0:

The pixel position on the image corresponding to the center pixel of the returned stamp. If either is None, they default to x0=floor(x+0.5) and y0=floor(y+0.5). (See above for why we don’t use round().) The peak* of the PSF on the returned stamp will be at (x-x0,y-y0) relative to the center pixel of the returned stamp.

“peak” assumes the PSF is radially symmetric. If it’s not, by “peak” read “center” or “fiducial point”.

Lots and lots of notes and examples to think through exactly what this means:

Algebra:

Define xc = floor(x + 0.5), yc = floor(y + 0.5). This is: the “closest integral pixel position” on the original image to where the PSF is being rendered. (It’s slightly different from the pixel position rounded to the nearest integer; see above.)
Define fx = x - xc, fy = y - yc ; both are in the range: [-0.5, 0.5).
Define midpix = stamp_size // 2: (so, for instance midpix=3 for a 7×7 stamp)
Given how we’ve defined the x and y parameters to this: function, on the original image, the peak of the PSF is at (x, y) = (xc + fx, yc + fy).
Pixel (midpix, midpix) on the stamp corresponds to (x0, y0): on the original image (given the definition of the parameters to this function).
If (x0, y0) = (xc, yc), then the “closest integral pixel: position” for the peak of the PSF on the stamp is (midpix, midpix).
In general, the “closest integral pixel position” for the: peak of the PSF on the stamp is (midpix + xc - x0, midpix + yc - y0). (If xc is 5 and x0 is 6, then the center of the stamp is to the right of the peak pixel on the stamp, so the peak pixel position is less than midpix.)
The peak position of the PSF on the stamp is: (midpix + xc - x0 + fx, midpix + yc - y0 + fy) (which is the same as (midpix + x - x0, midpix + y - y0)).
If we define (xrel=0, yrel=0) to be the peak of the PSF, then: (xrel, yrel) = (0, 0) is (midpix + xc + fx - x0, midpix + yc + fy - y0) on the stamp.
Therefore the center of the stamp, (midpix, midpix), is at: (xrel, yrel) = (x0 - xc - fx, y0 - yc - fy)
The center of the lower-left pixel of the stamp, (0, 0), is at: (xrel, yrel) = (x0 - xc - fx - midpix, y0 - yc -fy - midpix)

Examples:

For example: if you call psfobj.get_stamp(111., 113.), and if the PSF object as a stamp_size of 5, then you will get back an image that looks something like:

-----------
| | | | | |
-----------
| |.|o|.| |
-----------
| |o|O|o| |
-----------
| |.|o|.| |
-----------
| | | | | |
-----------

the PSF is centered on the center pixel of the stamp (i.e. 2,2), and that pixel should get placed on pixel (x,y)=(111,113) of the image for which you’re rendering a PSF. (Suppose you wanted to add this as an injected source to the image; in that case, you’d add the returned PSF stamp to image[111:116,109:114] (remembering that numpy arrays of astronomical images using all the defaults that we use in this software are indexed [y,x]).)

If you want an offset PSF, then you would use a different x0, y0. So, if you call psfobj.get_stamp(111., 113., x0=112, y0=114), you’d get back:

-----------
| | | | | |
-----------
| | | | | |
-----------
|.|o|.| | |
-----------
|o|O|o| | |
-----------
|.|o|.| | |
-----------

In this case, center pixel of the returned stamp corresponds to pixel (x,y)=(112,114) on the image, but the PSF is supposed to be centered at (x,y)=(111,113). So, the PSF is one pixel down and to the left of the center of the returned stamp. The peak of the PSF is at pixel (x-x0,y-y0)=(-1,-1) relative to the center of the stamp. If you wanted to add this as an injected source on to the image, you’d add the PSF stamp to image[112:117,110:116] (again, remembering that numpy arrays are indexed [y,x]).

If you call psfobj.get_stamp(111.5,113.5), then you’d get back something like:

-----------
| | | | | |
-----------
| |.|.| | |
-----------
|.|o|o|.| |
-----------
|.|o|o|.| |
-----------
| |.|.| | |
-----------

Because your pixel position ended in (0.5, 0.5), the PSF is centered on the corner of the pixel. The center of the stamp (x,y)=(2,2) corresponds to (floor(111.5+0.5), floor(113.5+0.5)) on the image, or (x,y)=(112,114).

If you call psfobj.get_stamp(111.5, 113.5, x0=111, y0=113) then you’d get back a stamp:

-----------
| | |.|.| |
-----------
| |.|o|o|.|
-----------
| |.|o|o|.|
-----------
| | |.|.| |
-----------
| | | | | |
-----------

Finally, to belabor the point, a couple of more examples. If you call psfobj.get_stamp(111.25, 113.0), you’d get back a stamp with the peak of the psf at (x,y)=(2.25,2.0) on the stamp image, with the center pixel corresponding to (x,y)=(floor(111.25+0.5), floor(113.+0.5)), or (111,113). You would add it to image[111:116,109:114], and the stamp would look like:

-----------
| | | | | |
-----------
| | |o|.| |
-----------
| |.|O|o|.|
-----------
| | |o|.| |
-----------
| | | | | |
-----------

If you call psfobj.get_stamp(111.25, 113.0, x0=110, y0=114), then you’d get a PSF back with the peak of the PSF on the stamp at (x,y)=(3.5,1.0), the center pixel corresponding to (x,y)=(110,114) on the image, and a stamp that looks like:

-----------
| | | | | |
-----------
| | | | | |
-----------
| | | |o|.|
-----------
| | |.|O|o|
-----------
| | | |o|.|
-----------

The peak of the PSF is at (x-x0,y-y0)=(1.25,-1.0) relative to the center of the returned stamp.

type x0:

int, default None

param y0:

The pixel position on the image corresponding to the center pixel of the returned stamp. If either is None, they default to x0=floor(x+0.5) and y0=floor(y+0.5). (See above for why we don’t use round().) The peak* of the PSF on the returned stamp will be at (x-x0,y-y0) relative to the center pixel of the returned stamp.

“peak” assumes the PSF is radially symmetric. If it’s not, by “peak” read “center” or “fiducial point”.

Lots and lots of notes and examples to think through exactly what this means:

Algebra:

Define xc = floor(x + 0.5), yc = floor(y + 0.5). This is: the “closest integral pixel position” on the original image to where the PSF is being rendered. (It’s slightly different from the pixel position rounded to the nearest integer; see above.)
Define fx = x - xc, fy = y - yc ; both are in the range: [-0.5, 0.5).
Define midpix = stamp_size // 2: (so, for instance midpix=3 for a 7×7 stamp)
Given how we’ve defined the x and y parameters to this: function, on the original image, the peak of the PSF is at (x, y) = (xc + fx, yc + fy).
Pixel (midpix, midpix) on the stamp corresponds to (x0, y0): on the original image (given the definition of the parameters to this function).
If (x0, y0) = (xc, yc), then the “closest integral pixel: position” for the peak of the PSF on the stamp is (midpix, midpix).
In general, the “closest integral pixel position” for the: peak of the PSF on the stamp is (midpix + xc - x0, midpix + yc - y0). (If xc is 5 and x0 is 6, then the center of the stamp is to the right of the peak pixel on the stamp, so the peak pixel position is less than midpix.)
The peak position of the PSF on the stamp is: (midpix + xc - x0 + fx, midpix + yc - y0 + fy) (which is the same as (midpix + x - x0, midpix + y - y0)).
If we define (xrel=0, yrel=0) to be the peak of the PSF, then: (xrel, yrel) = (0, 0) is (midpix + xc + fx - x0, midpix + yc + fy - y0) on the stamp.
Therefore the center of the stamp, (midpix, midpix), is at: (xrel, yrel) = (x0 - xc - fx, y0 - yc - fy)
The center of the lower-left pixel of the stamp, (0, 0), is at: (xrel, yrel) = (x0 - xc - fx - midpix, y0 - yc -fy - midpix)

Examples:

For example: if you call psfobj.get_stamp(111., 113.), and if the PSF object as a stamp_size of 5, then you will get back an image that looks something like:

-----------
| | | | | |
-----------
| |.|o|.| |
-----------
| |o|O|o| |
-----------
| |.|o|.| |
-----------
| | | | | |
-----------

the PSF is centered on the center pixel of the stamp (i.e. 2,2), and that pixel should get placed on pixel (x,y)=(111,113) of the image for which you’re rendering a PSF. (Suppose you wanted to add this as an injected source to the image; in that case, you’d add the returned PSF stamp to image[111:116,109:114] (remembering that numpy arrays of astronomical images using all the defaults that we use in this software are indexed [y,x]).)

If you want an offset PSF, then you would use a different x0, y0. So, if you call psfobj.get_stamp(111., 113., x0=112, y0=114), you’d get back:

-----------
| | | | | |
-----------
| | | | | |
-----------
|.|o|.| | |
-----------
|o|O|o| | |
-----------
|.|o|.| | |
-----------

In this case, center pixel of the returned stamp corresponds to pixel (x,y)=(112,114) on the image, but the PSF is supposed to be centered at (x,y)=(111,113). So, the PSF is one pixel down and to the left of the center of the returned stamp. The peak of the PSF is at pixel (x-x0,y-y0)=(-1,-1) relative to the center of the stamp. If you wanted to add this as an injected source on to the image, you’d add the PSF stamp to image[112:117,110:116] (again, remembering that numpy arrays are indexed [y,x]).

If you call psfobj.get_stamp(111.5,113.5), then you’d get back something like:

-----------
| | | | | |
-----------
| |.|.| | |
-----------
|.|o|o|.| |
-----------
|.|o|o|.| |
-----------
| |.|.| | |
-----------

Because your pixel position ended in (0.5, 0.5), the PSF is centered on the corner of the pixel. The center of the stamp (x,y)=(2,2) corresponds to (floor(111.5+0.5), floor(113.5+0.5)) on the image, or (x,y)=(112,114).

If you call psfobj.get_stamp(111.5, 113.5, x0=111, y0=113) then you’d get back a stamp:

-----------
| | |.|.| |
-----------
| |.|o|o|.|
-----------
| |.|o|o|.|
-----------
| | |.|.| |
-----------
| | | | | |
-----------

Finally, to belabor the point, a couple of more examples. If you call psfobj.get_stamp(111.25, 113.0), you’d get back a stamp with the peak of the psf at (x,y)=(2.25,2.0) on the stamp image, with the center pixel corresponding to (x,y)=(floor(111.25+0.5), floor(113.+0.5)), or (111,113). You would add it to image[111:116,109:114], and the stamp would look like:

-----------
| | | | | |
-----------
| | |o|.| |
-----------
| |.|O|o|.|
-----------
| | |o|.| |
-----------
| | | | | |
-----------

If you call psfobj.get_stamp(111.25, 113.0, x0=110, y0=114), then you’d get a PSF back with the peak of the PSF on the stamp at (x,y)=(3.5,1.0), the center pixel corresponding to (x,y)=(110,114) on the image, and a stamp that looks like:

-----------
| | | | | |
-----------
| | | | | |
-----------
| | | |o|.|
-----------
| | |.|O|o|
-----------
| | | |o|.|
-----------

The peak of the PSF is at (x-x0,y-y0)=(1.25,-1.0) relative to the center of the returned stamp.

type y0:

int, default None

param flux:

Ideally, the full flux of the PSF. If your stamp is big enough, and the PSF is centered, then this will be the sum of the returned stamp image. However, if some of the wings of the PSF are not captured by the boundaries of the PSF, then the sum of the returned stamp image will be less than this value.

type flux:

float, default 1.

rtype:

2d numpy array

Navigation

Related Topics

PSF¶

Note:¶

INDEXING IMAGES¶

PSF CENTERING FOR get_stamp¶