This looks great! We might want to think a little about how intelligent to make the detection algorithm / photometry functions -- e.g., do we want to detect and somehow notify the users about saturated or blended sources?

I don't see any discussion or capabilities for local background estimation - i.e. using an annular aperture to measure the local sky level and subtract that off of the counts within the source aperture. That's an essential capability for many use cases.

Will there be support for different image normalization conventions? In particular, it would be good to be able to handle the per-beam normalization used in radio and sub-mm/mm astronomy in addition to the usual optical convention.

It might be useful to add the normalization as metadata to NDData instances so that the code could take care of this automatically.

If you aren't familiar with this, in per-beam normalization the peak of the psf is the flux of the source (for unresolved objects), as opposed to the sum of all the pixels. This may seem strange (it did to me when I moved from optical to sub-mm), but it makes sense when you realize that pixels are an arbitrary user choice for radio/sub-mm/mm data.

@mperrin Local background subtraction could be done using:

r = 3.
ap = CircularAperture(r)
bkgap = CircularAnnulus(6., 8.)

total = aperture_photometry(im, (x, y), ap)
bkgavg = aperture_photometry(im, (x, y), bkgap, statistic='mean')

flux = total['flux'] - bkgavg['flux'] * np.pi * r ** 2

... and you'd also want to combine the errors in each. I agree there should be an easier way to do this, but the base functionality is there. perhaps a bkgap argument to aperture_photometry? That would work, but you'd have to allow the user to specify multiple statistics (one for the main aperture and one for the background aperture). Might get messy.

@kbarbary @mperrin - I think adding a background argument to aperture_photometry and psf_photometry would make sense, though I agree that we need to think about how the statistic would work. There's no reason we couldn't have statistic_background too.

Since you explicitly asked to not discuss implementation details at this point, I won't, hopefully the following question won't be too cryptic:

Wouldn't it make sense to move some of the underlying machinery needed to implement the photometry functionality in separate modules. Specifically:

• An astropy.photometry.Aperture could be a astropy.region.Region. A Region would e.g. also be useful to do a spatial selection on an unbinned event list.
• An astropy.photometry.PSF could be a astropy.???.Kernel. A Kernel would e.g. also be useful to convolve the image.

I guess applying a spatial filter to an event list or PSF-convolving an image is functionality we eventually want in astropy, but are (at least at the moment) not part of astropy.photometry, so it might make sense to be more modular from the start, to avoid backward-incompatible refactoring a year from now?

What about ds9 or CIAO region files?
Could you add support for reading / writing / exchanging apertures via strings or files?

@cdeil - thanks for your ideas!

FYI, I am waiting until Astropy v0.2 is ready to be released to start responding and implementing changes to the API document. Please feel free to continue leaving comments on this!

Just a note following a discussion with @joeharr4 - we need to support different interpolation techniques for the sub-pixel sampling. This can easily be done for the non-exact mode because we can use bilinear interpolation, or more generally scipy.ndimage.zoom which can easily oversample an image using spline interpolation of various orders. Whether we would be able to provide an 'exact mode' with interpolation remains to be seen.

For bilinear interpolation, interp2d would work.

@astrofrog I didn't exactly understand your note about interpolation. Is this in reference specifically to PSF photometry?

I was referring to aperture photometry - when we sub-sample, we currently simply subdivide the pixels into subpixels with the same value as the parent, so we are effectively choosing a finer grid and doing nearest neighbor interpolation. But you could also imagine defining the same fine grid, but instead using e.g. bilinear interpolation to find the values on the fine grid. Have a look at scipy.ndimage.zoom and try it out on an array to see what I mean. This exists in some other aperture photometry codes (I'll send you examples).

Ah I see. What we are currently doing is assuming the flux is distributed evenly within each pixel, whereas a bilinear interpolation would be guessing at how flux is distributed within a pixel based on the values of its neighbors. Is that right? (This may be a CCD-specific way of thinking about it)

Yes, that's correct!

It's important NOT to offer just any interpolation scheme. Bilinear is
flux-conserving, as is nearest-neighbor, but most/all others are not.
Obviously, losing flux eliminates the point of photometry.

I'm on vacation this week, but will have a closer look at this thread
next Monday or Tuesday.

--jh--

I've now got time to work on this again, so I've gone through the comments, and the main one that @cdeil and @perrygreenfield both raised is whether PSFs and apertures are defined in world or pixel coordinates. This is a fairly big issue, so I'm going to raise it on the list (and then we can get back to the more minor points)

I'd like to re-start the discussion on this photutils API spec.

@bsipocz will start her GSoC which is mainly about integrating photutils with astropy.coordinates, astropy.wcs, astropy.nddata and astropy.table two weeks from now. It's OK if we continue the photutils API discussion over the coming months as the implementation progresses and we realise what works well and what doesn't, but it can't hurt to start now.

I have the following questions:

• Is there a plan for a PixelCoordinates class or should the final photutils API use (x, y) tuples as input / output for pixel coordinates? IMO it's a bit inconsistent to use tuples for pixel and classes for world coordinates and the risk of users getting the photutils coordinate convention wrong and slightly shifted science results could be increased with (x, y) tuples. If others think a PixelCoordinates class is a good idea, does anyone have a reference to this in the coordinates or generalised WCS work that is going on or is this our responsibility to define it here?
• What do people think about storing the aperture position in the Aperture class? I would find this more natural than the current approach of having x, y separate from the Aparture class.
• @astrofrog has proposed to rename the annulus functions in Rename annulus functions photutils#34 ... please comment.

[newbie warning, possibly thinking orthogonally]
After advocating above for Aperture definition in world coords, I tried to describe photometry in its most general form, because I think it would be useful to think about it this way and then decide how specific to make these tools, rather than trying to back-engineer generality later.

Photometry: one desires the flux in some angular region of the sky, and range of wavelength. I assume this region has sharp boundaries (rectangular response function).

Real life: the sky brightness distribution is

1. convolved with the spatial instrumental response "beam"
2. convolved with the spectral response
3. convolved with the instrumental sampling function "pixel"
4. measured with some noise function - for photon counting detectors, that's the Poisson noise on a per-pixel measurement and one needs to know the gain (+readnoise etc), but that's not always the case - for a synthesis image its wave noise / system temperature on a per-beam measurement, so in the general case for noise calculation, the noise can't be assumed uncorrelated per pixel. In general I think one needs a noise map as a function of position, with units (counts/pixel for a CCD, Jy/bm for a synthesis image), and then the photometry routine calculates how many independent measurements there are in the Aperture and propagates uncertainty accordingly.

[re] 3. I believe that only the case of the space-filling, contiguous rectangular response function pixel has been considered so far here - in effect the code to deal with fractional pixels is just deconvolving the pixel response from the image, isn't it?
I raise for consideration whether the assumption of that specific pixel response function may be overly restrictive. Although I've never done this myself, I'm told that accurate CCD photometry does account for the non-rectangular pixel response. Another case is synthesis imaging, in which the pixel sampling function is a delta function. If a base class can be straightforwardly designed with this deconvolution in mind, then inherited classes could deal with different deconvolutions?

[re] 2. I would have thought complex spectral response to be the most "special-case" part of this, and that we could design tools assuming effective delta-function or "equivalent wavelength" spectral response, but @cdeil has already offered a more complex use case.

[re] 1. As far as I can tell this hasn't been considered yet (i.e. no aperture correction) except implicitly in the normalization of the PSF for PSF photometry.

So... In my mind, there's an Aperture in world space (3d), and then

• an Image which ideally knows WCS, beam and pixel response functions, measured intensity and noise at discrete locations in pixel space
• or an EventList with position and response functions..

And although Tom doesn't want implementation to bleed into the discussion, I guess the Image would know, when presented with an Aperture, how to create a selection function a.k.a. mask/weight image in pixel coordinates, and apply it with unit conversion?

Maybe this is too far abstract and afield, in which case I'll take silence for an answer :)

@indebetouw Thanks for joining this discussion! I'll have to think about your last comment some more and try to figure out how it should be reflected in this API document.

Everyone: there has been lots of discussions, but the actual document hasn't been updated for over a year. What would be most helpful at this point IMO would be:

• Pull requests against this pull request by @astrofrog with concrete changes / additions to be made.
• Code examples (e.g. in gists or IPython notebooks, doesn't have to work) demonstrating your use case that is not well supported by the current API.

I'll try and find the time to do this myself for my use case (gamma-ray event data with RA, DEC, Energy detected for each photon) soon.

I've been lurking on this discussion since the outset, and I've been
debating on whether to say something or just give up.

Here's the basic gist of what I'm going to say:

Keep the default case simple and free of assumptions about the use case.
Layer routines that are good for specific purposes on top of that,
ideally without repeating code.

From the discussion, the photometry system you folks are putting
together seems like it might be useful for some set of idealized
observations that you as individuals happen to do, and it might appeal
to CS people because of its cool and flexible design, but if it's much
more complicated than specifying a center in pixel coordinates and an
aperture and two annular radii, only a limited number of people will use
it. And, that specification has to be available in a simple, procedural
way as arguments to a function call. And, it cannot depend on anything
in the FITS header. And, it shouldn't do anything extra or unexpected.

This isn't to say that you can't offer WCS or that an OO way of calling
it shouldn't be available; of course they should. But, if the most
basic way of calling it is much more than:

(flux, uncert [, raw, skyave, npix, napbad, nskybad]) =
apphot(image, (y, x), aprad, skyin, skyout,
maskim=False, uncertim=False, gain=False, fulloutput=False)

where [] indicates the optional stuff that fulloutput=True enables, then
it's too complex and people will just continue to use their own homebrew
routines. And if it's much less, it makes too many assumptions and
again, people will stick with what they have, or even write simple
apphot routines that fit their purposes in order to avoid yours.

So, I'd take the Matplotlib approach and make sure that such a routine
with pretty much these args exists, and that the internal design isn't
so elaborate that it makes this difficult to implement.

Here are specific concerns:

It has to work on plain 2D numpy arrays without a FITS header and
without going through an extra setup routine to make some kind of
special image object. The images may come from a non-FITS source, such
as a simulation. The dataset could be a 3D image cube with one
representative FITS header (e.g., Spitzer subarray sets), and the WCS
info may be wrong for all but one image. If you have 1e4 or 1e5 images
in a "cube", as time-series datasets often do, you are probably not
carrying the headers around, as they may be bigger than the images.
Likely they weren't even generated by the instrument.

I encounter very few people who use WCS in any way. This isn't to say
that WCS isn't useful to many people, but it does mean that there are
TONS of astronomers (pretty much everyone in planetary science, and
much/most of the exoplanet characterization community) who find no
benefit in WCS. Either our objects move, or the precision of the
nominal WCS solution in an image is too poor to do the job at the level
we need it done (0.02 pix centering accuracy, in the case of
high-precision Spitzer/HST work).

Likewise, centering is a separate issue from photometry. How best to do
it depends on the application, resolution, and S/N. So, keep centering
in a separate set of routines. Obviously, the output of those routines
needs to be in the same format as the input of photometry.

You could have a center=False option that could be reset by the user to
the name of a centering routine, of course, but that should not be the
default. Likewise with WCS, a more OO way of accessing the routines, or
arrays of coordinates or 3D stacks of images.

Keep the default case simple and free of assumptions about the use case.
Layer routines that are good for specific purposes on top of that,
ideally without repeating code.

Thanks,

--jh--
Prof. Joseph Harrington
Planetary Sciences Group
Department of Physics
PSB 441
4000 Central Florida Blvd.
University of Central Florida
Orlando, FL 32816-2385
jh@physics.ucf.edu
planets.ucf.edu
skype: joeharr4

+1 to everything Joe said.

Although I've been listening in, I haven't been active in photutils discussions lately, because I've mostly taken Joe's option of "just give up."

Tangentially, I have a lot I could say about performance, but a summary would be:

• It's very important to some people (like me)
• a design that Joe is advocating would probably be easier to optimize
• Since writing it, I've realized that even the existing code in photutils could be much faster if it was designed differently.

I'd be happy to talk more about performance. On the other hand, if no one else cares about it, I won't say any more.

I'll have a think about these comments, but I know @cdeil is already planning to open a PR to simplify some of the items in this API proposal. I'm surprised about the WCS comments - surely people would use the ability to specify source positions in e.g. equatorial coordinates?

I think the general tool needs specification in both world or pix coords. The question in my mind is whether one has a low-level tool that's not wcs-aware, and then something above which does the wcs transformations and calls the lower level code, or whether one builds the Aperture to handle either case and transform between them.
If one has a lot of data from one instrument (and especially if that instrument is bad at writing wcs headers), then one doesn't want to do the wcs x-forms multiple times, if at all.
If one has a multi-instrument dataset, one likely has to work in world coords for consistency.

@indebetouw - I agree, I think that at the very core, there should of course there should be a low-level function that is not WCS-aware. I don't think that providing the optional convenience of passing world coordinates to the aperture class will affect performance. @cdeil and I will be updating the document in the next few days.

I'll make detailed comments after the update @astrofrog mentioned, but one API suggestion for the short term: maybe lead with a few high-level functions that tie together the underlying implementation details. Some of the concerns @joeharr4 expressed are explicitly addressed already in the API but pretty far down -- there is no intent as far as I can tell to restrict the input data to NDData or CCDData or some other Image class.

My experience with ccdproc was that I wanted the initial API as detailed as possible to make for a clear initial coding goal, but most users won't need to know about most of the API.

Aperture with annulus: I do think that there needs to be a convenience function or two to make the common case of using an annulus to determine the background easy. Getting the common case of many apertures and many stars to work with the current implementation is a little tricky; it requires either looping or diving into the rules of numpy array broadcasting. It ends up being a small number of lines of code, though.

WCS/pixel space: Both/either please :). Most of the time I'll prefer to give a position in world coordinates, but will have instances where I need to work with pixels instead.

It would be really nice if this does end up supporting 3D aperture spectrophotometry as Christoph & others suggest (I think) :-). I believe I still haven't run into a good tool for doing that yet. Of course that might be for a later version.

@joeharr4 I liked your suggestion so much, I've implemented something pretty close to it in my SEP project. See https://github.com/kbarbary/sep-python under "aperture photometry" under "Usage". Beta testers wanted.

The idea is that SEP will hopefully be used in photutils for the "backend".

@astrofrog - Given that most of this is now implemented in photutils 0.1, I'm not sure it's useful to finish this document. Close this issue and add a link to it to the photutils docs so that the initial API spec and this discussion is preserved in case someone proposes concrete API changes in the future?

@cdeil - maybe I can add a note at the top that says that is is non-binding and then merge? Otherwise it will be lost if I delete my fork in future.

+1 to adding a note at the top saying this is somewhat outdated and non-binding and merging it.
This makes it clear that further discussion / work should go into concrete photutils issues / PRs.