RomanDatamodelPSF

class snappl.psf.RomanDatamodelPSF(_parent_class=False, stamp_size=None, **kwargs)

Bases: PSF

A PSF pulled from the gridded oversampled PSF stored in images from roman_datamodel.

TODO WORRY : figure out how these things are centered! The example I’m looking at has a 361×361 image for the PSF, with an oversample factor of 4. Turns out that 361 is not an integer multiple of 4. So how is this really gridded???? I don’t know. Current code assumes that the oversampling is such that the center of the image it gets is the center of the oversampled PSF, so the center pixel of an aligned image-resolutin PSF will have its borders drawn through the edge of the oversampled pixels, not on pixel edges.

(Anecdote: looking at the results of one get_stamp, it “seems OK”.)

TODO: Allow selection of the SED star model used. Right now it’s just taking the second one (index 1, which is G2V… at least, it was in the example I looked at). This has implications for the interface of the base PSF class, however, so more work is required.

WRITE DOCSTRING

Parameters:

stamp_size (the stamp_size of the stamp that will be returned by get_stamp. Must be odd.)

Attributes Summary

stamp_size

The size of the one side of a PSF image stamp at image resolution.

Methods Summary

getPhotutilsPSF([x, y])	Return a photutils.psf.SOMETHING that can be the psf_model parameter of photutils.psf.PSFPhotometry
get_stamp([x, y, x0, y0, flux])	Return a 2d numpy image of the PSF at the image resolution.

Attributes Documentation

stamp_size

Methods Documentation

getPhotutilsPSF(x=None, y=None)

Return a photutils.psf.SOMETHING that can be the psf_model parameter of photutils.psf.PSFPhotometry

You really only want to use the result of this for feeding to photutils. If what you want is an image of the PSF, use get_stamp() instead.

Unless the subclass implementes something else, this will be a photutils.psf.ImagePSF, and will be just for one position. If you don’t give x and y here, it will use the x and y that were given to construct the PSF object.

Some subclasses will ignore x, y and return a Photutils PSF object that can be used at any x, y. Read the subclass documentation to find out what is what. If you’re writing code to work with all subclasses, and you are dealing with a spatially variable PSF, then be careful with what you do with the return here.

NOTE AND WORRY : photutils uses the term PSF to mean both PSF and PRF in different places. Just like everyboyd else. THIS is returning a PRF (if the code is right), that is, a thing that has been convolved with a tophat pixel (or, ideally, something even more sophisticated, but no promises). However, we call it a PSF, because it’s the thing you give to photutils PSF photometry. But photutils has started to try to distinguish, so their documentaiton is going to be scary and inconsistent and we’ll have to think.

Subclasses may want to override this.

Parameters:

• x (float, default None) – Position on the original image to render the PSF for. If not given, uses the x and y that the PSF object was initialized with. Because photutils expects PSFs centered on a pixel, x and y should have no fractional part. (The default implementation explicitly rounds (sort of…) to the nearest whole number for both x and y; properly written subclasses will do the same thing.

• y (float, default None) – Position on the original image to render the PSF for. If not given, uses the x and y that the PSF object was initialized with. Because photutils expects PSFs centered on a pixel, x and y should have no fractional part. (The default implementation explicitly rounds (sort of…) to the nearest whole number for both x and y; properly written subclasses will do the same thing.

• imagesampled (bool, default True) – This one is a bit complicated. It’s ignored by some subclasses. The default implementation is to return an ImagePSF sampled at the image scale, and to ignore imagesample. Some subclasses (e.g., OversampledImagePSF) will return an oversampled photutils.psf.ImagePSF if you set imagesampled to false. Other subclasses return other things for which image scale may or may not be relevant. And, it’s possible that other subclasse will return something oversampled even if you say imagesampled=True. In any event, because what you get back is intended to be fed to photutils, you probably shouldn’t think about it too hard. (TODO: we might want to set imagesampled to False so that classes don’t by default make something that will make photutils behavior worse! This may cause a bunch of code to break, maybe quietly, because the assumed default changed.)

get_stamp(x=None, y=None, x0=None, y0=None, flux=1.0)

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