Proposed Photometry API

[astrofrog]

As discussed at the Astropy coordination meeting at STScI (https://github.com/astropy/astropy/wiki/Minutes-for-AstroPy-Meeting-2012), having the ability to do photometry within Astropy would be a great feature.

The attached document presents a proposed API for photometry in Astropy. Photometry functionality is currently being developed in photutils (http://github.com/astropy/photutils) but we intend that this would eventually be merged in to astropy.photometry.

At this time, we solicit comments from the Astropy developer community about this proposed user API. Please let us know if any crucial functionality is missing, or if you think the API could be improved or clarified.

Please note: unless they have a direct impact on the user API, implementation details should not be discussed at this stage.

This document was prepared by @astrofrog, @kbarbary, @bretonr, and @stargaser, and has already been through several iterations before the current version.

[adrn]

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?

[mperrin]

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.

[mperrin]

From the way things are written here, it sounds like you cannot specify the gain if you're using a NDData object? I haven't followed the development of NDData closely enough to know, is there an official, unambiguous way that gain will be specified for all NDData objects? Otherwise you'll want to still be able to set that parameter here, I think.

[kbarbary]

• I think the photometry routines should accept only Image objects (Image would be a subclass of NDData). Image objects would have a gain property, which would default to 1.
• If an Image object is passed, it would still be useful to allow the user to specify these additional parameters. They would override the properties of the Image object if specified.

Note that in the current implementation, gain can be a scalar OR an array (matching dimensions of the data).

[astrofrog]

Reactions were pretty negative to my suggestion of an Image sub-class on the list a little while back (https://groups.google.com/d/topic/astropy-dev/ReCISxQJgQQ/discussion) but I guess this will be a good time to revive this.

[kbarbary]

I forgot about that discussion! There were good points on both sides... I didn't mean to restart that discussion here. Regardless, we can define a standard place in NDData or Image, where the gain should be defined. Either in a gain property or in the metadata.

[astrofrog]

gain is something specific to CCDs, so it could belong in a CCDImage class.

[indebetouw]

Looking forward to proper uncertainty calculations, should there be a way to distinguish between unknown gain and gain=1. e.g. default gain=None, assumed=1 with a warning when used in calculations?

[cdeil]

In the meantime ccdproc has emerged which has a CCDData class as an NDData sub-class ... I'm not sure if photutils should be aware of CCDData contents or if it should be more of a toolbox that e.g. sums arrays over apertures that is called by ccdproc and other packages.

I think I'd prefer the "toolbox" approach so that photutils becomes usable e.g. also for count data in gammapy, but for very common cases like CCD data it could also be OK if photutils is CCDData-aware?

@mwcraig @crawfordsm What do you think how ccdproc and photutils should work?

[mwcraig]

I'd like ccdproc and photutils to work well with each other, of course, but I think we have to assume there will be users who want or need to use one but not the other.

CCDData currently does not store a gain; the only things that distinguish it from plain NDData right now are that the uncertainty is a StdDevUncertainty and it requires that a unit be set. The gain isn't even required to be in the metadata (and is currently ignored even if it is in the metadata, actually); when we need it, there is a function argument for providing it.

[crawfordsm]

Well, the idea of inheriting from NDData was so that we could pass around the objects from one tool to the next, so I would hope that the CCDData object would be usable by this class. As @mwcraig mentions, there shouldn't be anything in CCDData that should prevent it from being used vs. any other NDData object. So I might change the above in the API to something like 'if an object is inherited from an NDData, it can be passed to the function. ' rather than just an NDData object.

[crawfordsm]

I'd also add that if the photometry tools haven't actually been written yet, it might be worthwhile to re-visit the idea of having an ImageData class. Although we have gone back and forth on this debate, having a CCDData class has made it much easier to write the code because: a) we didn't have to specify parameters for every function, b) any checking on what is needed for the object can be done in the CCDData class instead of repeating in every function, and c) the correct object with the needed properties was always being passed to the functions.

Of course, we are the only ones to use ccdproc so far, so although it might be easier to code, people may hate using it in that format, but I thought I'd share my experiences so far with ccdproc.

[aconley]

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.

[perrygreenfield]

What frame are the apertures defined in? Pixel coordinates or world coordinates? There are telescope/detector combinations that this makes a significant difference. Does specifying pixels or world coordinates select the frame the aperture is computed for?

[astrofrog]

This is a very good point that will matter for any non-circular aperture/PSF. How do you think it should work?

[perrygreenfield]

It may even matter for circular apertures (e.g., for ACS, the plate scale is on the order of 10% different in two different directions (iirc). A circle in one frame is an ellipse in the the other.

Offhand, I would think physically, one wants fixed apertures in world coordinates. E.g., that's where we want the constant angular area.

[cdeil]

I'm not sure, but I think there's also "sky coordinates" which would be different from "pixel coordinates" and "world coordinates". E.g. imagine you have a survey map like this one from Fermi and you want to do aperture photometry with circles of 1 or 10 deg radius on the sky or PSF photometry with Gaussians of width 1 or 10 deg on the sky.

So the main use case I have in mind with "sky coordinates" is that you have an aperture or PSF which only depends on the spherical distance on the sky.

I think computationally photometry with "sky coordinates" would be similar to photometry with "world coordinates", because for each given position you first have to compute the binary aperture mask or weights image, whereas for "pixel coordinates" you just shift one given aperture mask or PSF weights image around the image.

Am I making any sense?

I'm not saying astropy.photometry should support "sky coordinates" necessarily, I just wanted to mention it as one possible use case and that the API should be clear here what exactly it can an cannot do.

[astrofrog]

@cdeil - could you clarify what you see as the difference between sky and world coordinates? In WCS-speak, sky coordinates are just a subset of world coordinates (i.e. world coordinates that are positions on the sky), but maybe you are thinking of something else?

[cdeil]

They might or might not be the same. The problem is that I don't know what a CircularAperture(20 * u.degree) means say for the position (lon, lat) = (0, 60). In sky coordinates I mean all pixels that have spherical distance less than 20 deg from that coordinate. Is that the same meaning for world coordinates?

http://www.astro.rug.nl/software/kapteyn/plot_directive/EXAMPLES/mu_skypolygons_allsky.hires.png

[astrofrog]

My thinking was that it is the polygon that is defined as being 20 degrees on the sky from the position specified (as in, spherical distance). Just to make sure I understand, what other interpretation do you think this could have?

[cdeil]

For your example EllipticalAnnulusAperture(1.5, 2.0, 1.0, np.pi / 3.) above and the "world" or "sky" coordinate case, how is the orientation on the sky actually defined? I guess it depends on the coordinate system used in my image, e.g. if it is Galactic or Equatorial, right?

Does the orientation / shape of that aperture also depend on the FITS projection of my image, e.g. AIT or CAR?
If the answer here is no, then I agree, there is no third sky case besides pixel and world.

[astrofrog]

@cdeil - I want to revive this discussion. To answer your question about whether the shape of the aperture depends on the FITS projection of the image - I guess it could, since the world-to-pixel coordinates change. I guess I still don't understand the distinction between sky and world though. The aperture is just defined as a polygon on the sky with the desired properties, and mapped onto the pixel space. So say for the case of a large circular aperture, this aperture is defined by the set of points at a given spherical distance from the position considered. I then map this into pixel coordinates. This will depend on the details of the projection, pixel orientation, etc. So would you call this sky or world? And what is the other case?

[indebetouw]

For me its important to have apertures defined in world coords. The conversion to pixel coords is then "trivial" with wcslib, if one allows an arbitrarily shaped aperture. I've already semi-volunteered to write PolygonAperture, which I think as long as one chooses enough vertices will cover this use case.
Overall, defining apertures in world space is potentially safer/better, since one could imagine some image with distortion, where the scientifically correct thing to do for objects throughout the field would be the same aperture size on the sky, not in number of pixels.
Defining apertures in world space also makes performing matched photometry on different images natural instead of complicated.

In my mind, the aperture is a world-coord entity, and then that entity plus a specific image wcs yields an aperture in pixel space. The pixel coords of an aperture only have meaning in the context of that specific image, whereas the aperture itself is more universal, like a Platonic ideal.

[kbarbary]

@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.

[astrofrog]

@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.

[cdeil]

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?

[cdeil]

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

[astrofrog]

@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!

[astrofrog]

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.

[astrofrog]

For bilinear interpolation, interp2d would work.

[kbarbary]

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

[astrofrog]

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).

[kbarbary]

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)

[astrofrog]

Yes, that's correct!

[joeharr4]

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--

[astrofrog]

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)

[cdeil]

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.

[indebetouw]

[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 :)

[cdeil]

@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.

[joeharr4]

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

[kbarbary]

+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.

[astrofrog]

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?

[indebetouw]

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.

[astrofrog]

@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.

[mwcraig]

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.

[jehturner]

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.

[kbarbary]

@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".

[cdeil]

@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?

[astrofrog]

@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.

[cdeil]

+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.