Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Measurements Upgrades: coordinates, bitflag, psf photometry #288

Merged
merged 16 commits into from
May 30, 2024

Conversation

guynir42
Copy link
Member

@guynir42 guynir42 commented May 23, 2024

I've neglected to add a few things that should have been properly propagated from cutouts to measurements or calculated for the new measurements:

  • the coordinates should include the centroid offset for each measurement.
  • the badness bitflag that gets propagated from Image all the way down to Cutouts needs to also be propagated to measurements (and I've added a disqualifier cut that doesn't save Measurements that have a non zero bad flag).
  • PSF photometry: I've added this part using the new image's PSF (rather than the ZOGY PSF which isn't able to get sub-pixel shifted clips).
  • Forced photometry on cutouts: each measurement is able to give the photometry (in any aperture, or PSF) for any ra/dec that is close to the position of the cutouts. This doesn't replace needing to do forced photometry on images where there was no detection and thus no cutout.

@guynir42 guynir42 changed the title Measurements Upgrades: coordinates and bitflag Measurements Upgrades: coordinates, bitflag, psf photometry May 28, 2024
@guynir42 guynir42 marked this pull request as ready for review May 28, 2024 12:08
@guynir42 guynir42 requested review from rknop and whohensee May 28, 2024 12:08
Copy link
Contributor

@rknop rknop left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some questions and comments.

(self.datasize - self.imsize) // 2 : (self.datasize + self.imsize) // 2,
(self.datasize - self.imsize) // 2 : (self.datasize + self.imsize) // 2,
]

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This does raise the issue of doing photometry in a radius that's larger than the image size. This should probably return an error, or at lest a warning.

(When did this come up?)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've been using fairly small cutouts (25 pixels) and some of the default apertures are really big. I think we talked about changing them.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm. This is pretty serious downside on doing photometry on cutouts rather than directly on the sub image. To really do it may require cutouts to be bigger than we would otherwise want them to be.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't really see why we need to have such large apertures.

Default is False.
For each aperture, will measure and reposition the centroid
this many times before moving on to the next measurement.
Default is 2.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For clarity: does it do iterations measurements, updating the position after each one, and then do a final pass, for a total of iterations+1 measurements? Or does it do iterations-1 measurements, updating the posititon after each one, and then make a final measurement, for a total of iterations measurements?

Copy link
Member Author

@guynir42 guynir42 May 28, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It will do number of radii * iterations + a final run of number of radii

if denominator == 0:
denominator = 1.0
elif abs(denominator) < 1.0:
denominator = 1.0 * np.sign(denominator) # prevent division by zero and other rare cases
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why only when abs(denominator) < 1.0?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I chose an arbitrary minimum value but it can be made much smaller if you think we might see such values.

photometry['moment_xx'] = cxx
photometry['moment_yy'] = cyy
photometry['moment_xy'] = cxy
# for each radius, do 1-3 rounds of repositioning the centroid
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm confused about this.

It does this recentering for every radius, but unless I'm reading something wrong, it throws out the centroids for everything except the last radius. There's just a single cx variable, so it's going to get overwritten every time thorugh the loop over radii; the saving of cx and cy to best_cx and best_cy happens after the radius loop.

Why bother doing this for all of the radii if the results aren't going to be kept for anything but the last radius? (I recognize that it's not identical, because the results of previous radii are used as starting points for later radii, but I'm guessing that unless iterations=1, that's going to have a very small effect.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The goal is to slowly lock in the correct position.
When the aperture is big, it has a good chance of having the PSF inside it's area, but it also has a lot of noise that might make it less accurate. Then a little bit smaller aperture has a good chance of still having the star but will improve the centroid position. The final (smallest) aperture gives the best centroid and we use that for actually measuring the flux and second moments.

# get the subtraction PSF or (if unavailable) the new image PSF
psf = self.cutouts.sources.image.get_psf()
psf_clip = psf.get_clip(x=image_pixel_x, y=image_pixel_y)
offset_ix = int(np.round(offset_x))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps something to worry about -- in models.psf, PSF.get_clip() has the following code:

        xc = int( np.floor(x + 0.5) )
        yc = int( np.floor(y + 0.5) )

Usually this should be the same as int(np.round()), but it won't when offset is exactly +0.5. (Or, for numbers bigger than ±1, when ( ( offset is even ) and ( offset > 0 ) and ( frac(offset)==0.5) ) or ( ( offset is odd ) and (offset < 0 ) and (frac(offset)==0.5) ).)

We should probably just edit PSF.get_clip() to use np.round().

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK I removed the 0.5 from get_clip.

The sigma to use for the clipping of the measurements. Default is 3.0.
iterations: int, optional
The number of iterations to use for the clipping of the measurements. Default is 3.
measurement_list_kwargs: dict, optional
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thinking about this, and looking at get_measurements_list(), I'm thinking that we should probably require a prov_hash argument to get_measurements_list() (and thus here as well).

You have a default, which is as sane as any default could be. (At least, until we have some sort of version tagging table where we can tag provenances with names; at that point, we could have a 'default' tag, and then the provenances tagged that way are what would be used as default.) What I'm worried about, though, is somebody in the future looking at this code and editing it, and seeing the methods that are available, and assuming that the defaults are going to be the right thing to use unless you know otherwise. Provenances are a complicated thing, though, and it's very easy to come up with situations where the default you have isn't the right one. Given that provenances are complicated, if we make the API clear that you have to understand them well enough to do the right thing, it will force people to think about it rather than making a call that is legal as defined that will return something different from what they were expecting. (And, in such a way that they wouldn't know it was different from what they were expecting.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I generally agree: this is a place that might surprise the user (which is bad in general). But I am also not happy to have provenance hashes as something you MUST know to use this function.
Right now it does something reasonable and can let you use it without getting into the habit of choosing provenances and copying hashes.
I don't know, I can go either way on this.

points = SkyCoord(ra, dec, unit='deg')

ra_mean = np.nansum(ra * flux) / np.nansum(flux)
dec_mean = np.nansum(dec * flux) / np.nansum(flux)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about the case where only one of (ra,dec) is NaN, and flux is not NaN? I suspect this will be rare (it may not even ever happen). But, in case it does, might be better to filter so that ra_mean doesn't include measurements where dec is NaN (and likewise for dec).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

...In fact, it's bad if flux is NaN but ra is not NaN, or vice versa. Then the normalization is wrong.

Probably need something like

good = ( ~( np.isnan( ra ) | np.isnan( dec ) | np.isnan( flux ) )
ra_mean = ra[good].sum() / flux[good].sum
dec_mean = dec[good].sum() / flux[good].sum

(Might also want to put in a flux > 0 filter, or maybe flux/dflux > 5..)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah that works. I'm going with flux/fluxerr > 3.0, unless you feel it really should be 5.0

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

3.0 is fine

models/objects.py Outdated Show resolved Hide resolved
models/objects.py Show resolved Hide resolved
flux_small_aperture = m.get_flux_at_point(m.ra, m.dec, aperture=1)
flux_large_aperture = m.get_flux_at_point(m.ra, m.dec, aperture=len(m.aper_radii) - 1)
flux_psf = m.get_flux_at_point(m.ra, m.dec, aperture=-1)
assert flux_small_aperture[0] == pytest.approx(m.flux_apertures[1], abs=0.01)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These verify that a new run of get_flux_at_point returns what is expected.

However, looking at the fixtures, it looks like the values inside ptf_datastore are created using a Measurer, which (I believe-- I haven't dug through the code to verify) will have itself been done through a get_flux_at_point call.

Are there any explicit tests of get_flux_at_point itself? (That is, make sure it's doing what we expect, not that it does the same thing every time.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The PSF flux is calculated using the same code, but the aperture flux is gotten through the photometry module.
I have some tests for that code in test_measuring and it seems to be doing the right thing, but maybe I can add something to test the PSF photometry over there.

models/cutouts.py Outdated Show resolved Hide resolved
Copy link
Collaborator

@whohensee whohensee left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nothing jumped out at me besides a few questions for my own sqlalchemy/general python knowledge. I'm not entirely familiar with the exact details of the math in all of these, but since I've been lookin at these objects for the past few days I was able to mostly follow. Generally looks good to me.

@guynir42 guynir42 requested a review from rknop May 30, 2024 12:22
@rknop rknop merged commit efab173 into c3-time-domain:main May 30, 2024
6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants