# Write down astrometry model

XMLWordPrintable

#### Details

• Type: Story
• Status: Done
• Resolution: Done
• Fix Version/s: None
• Component/s: None
• Labels:
• Story Points:
8
• Epic Link:
• Sprint:
Per-CCD Sprint 1
• Team:
Alert Production

#### Description

Just write down the model carefully, identifying free parameters, labeled cameraGeom systems, and FITS WCS frames like "intermediate world coordinates". Then take derivatives w.r.t. our desired parameters. Is the model linear? Can we make it linear via some transformation, or by fitting for sky-to-pixels instead of pixels-to-sky first? This may be useful for DM-19734 and DM-19735 as well.

#### Activity

Hide
Jim Bosch added a comment -

I think that we can replace the code using the metadata in cameraMapper.py with calling makeSkyWcs (DM-19735 version) using the VisitInfo boresight and rotator angle (possibly transforming rotator angle to be an appropriate orientation, but it may be correct already).

I'm generally very much in favor of using that makeSkyWcs overall instead of makeDistortedTanWcs for initial WCSs.; it's just much easier to understand what's going on  But we may want to at least have an option to do something in between the two current approaches, in which we use the CRVAL and CD matrix along with cameraGeom to infer the boresight location and rotation angle, in case some cameras don't provide that accurately directly in their headers.  On the other hand, in that case, maybe it'd be astro_metadata_translator's job to do that in order to correctly populate VisitInfo?  (Tim Jenness may have an opinion on this, or on how reliable we can generally expect header boresight and rotation to be).

If we do this, we probably should recompute the crpix/crval to recenter them on the chip instead of being relative to the boresight, though it may not matter in practice.

I don't think you want to do this - if you're using the cameraGeom distortion and it's good enough for a decent WCS (which I think we have to assume if we go this route), then there's nothing to be gained, and you're complicating and code that wants to use the internal coordinate systems in those WCSs significantly by breaking any relationship that IWC (which is defined by CRVAL/CRPIX) and FIELD_ANGLE.

Show
Jim Bosch added a comment - I think that we can replace the code using the metadata in cameraMapper.py with calling makeSkyWcs ( DM-19735 version) using the VisitInfo boresight and rotator angle (possibly transforming rotator angle to be an appropriate orientation , but it may be correct already). I'm generally very much in favor of using that makeSkyWcs overall instead of makeDistortedTanWcs for initial WCSs.; it's just much easier to understand what's going on  But we may want to at least have an option to do something in between the two current approaches, in which we use the CRVAL and CD matrix along with cameraGeom to infer the boresight location and rotation angle, in case some cameras don't provide that accurately directly in their headers.  On the other hand, in that case, maybe it'd be astro_metadata_translator's job to do that in order to correctly populate VisitInfo?  ( Tim Jenness may have an opinion on this, or on how reliable we can generally expect header boresight and rotation to be).   If we do this, we probably should recompute the crpix/crval to recenter them on the chip instead of being relative to the boresight, though it may not matter in practice. I don't think you want to do this - if you're using the cameraGeom distortion and it's good enough for a decent WCS (which I think we have to assume if we go this route), then there's nothing to be gained, and you're complicating and code that wants to use the internal coordinate systems in those WCSs significantly by breaking any relationship that IWC (which is defined by CRVAL/CRPIX) and FIELD_ANGLE.
Hide
Tim Jenness added a comment -

Regarding accuracy, it really depends on the telescope. Usually boresight rotation is pretty accurate because it's reading out directly from some hardware.

Regarding re-centering, I really don't think we should do that either. If you move CRPIX to the center of the chip you have to translate the polynomial distortion model to the new origin and that can introduce problems. What do you gain by doing that?

Show
Tim Jenness added a comment - Regarding accuracy, it really depends on the telescope. Usually boresight rotation is pretty accurate because it's reading out directly from some hardware. Regarding re-centering, I really don't think we should do that either. If you move CRPIX to the center of the chip you have to translate the polynomial distortion model to the new origin and that can introduce problems. What do you gain by doing that?
Hide
John Parejko added a comment -

What do you gain by doing that?

A concern we had was that the projection from IWC (tangent plane) to sky (spherical) could cause problematic distortions for the edge chips, if crpix/crval was defined at the boresight.

We will have to modify crval though, to correct for mis-centering, so the fitter will have to be smart about converting the matched refcat into crval at the boresight.

Show
John Parejko added a comment - What do you gain by doing that? A concern we had was that the projection from IWC (tangent plane) to sky (spherical) could cause problematic distortions for the edge chips, if crpix/crval was defined at the boresight. We will have to modify crval though, to correct for mis-centering, so the fitter will have to be smart about converting the matched refcat into crval at the boresight.
Hide
John Parejko added a comment -

# Design/change sketch for AstrometryTask updates

I think that the below encompasses all of the necessary work to implement the new single frame astrometry design. John Swinbank: is this sufficient to close out this ticket? I'm filing a new ticket to start the work on the first part.

## Reading raw Exposure changes

• This work is going to go into the Gen2 butler world. Converting it to gen3 butler will be a separate job, but could be made easier by keeping all of the relevant changes in free functions that we can just use directly.
• The initial SkyWcs is assigned in exposureFromImage (called from _standardizeExposure, called from std_raw, so this all happens when the raw is loaded by the butler via std_raw).
• exposureFromImage first calls makeSkyWcs(metadata), and only later builds the VisitInfo. We will have to reorder that so that VisitInfo gets built first, then makeskyWcs(pixelsToFieldAngle, orientation, flipX, boresight) is called later, using the camerageGeom to specify pixelsToFieldAngle and flipX, and the VisitInfo to specify orientation and boresight.
• Need to write a separate function to call from exposureFromImage (createInitialSkyWcs(visitInfo, detector)?) that looks similar to lsstCamMapper.getWcsFromDetector.
• We may have a problem due to when metadata gets stripped. Jim mentioned this gen3 commit as an example. However, if we don't actually need the FITS WCS metadata at all (which we won't with the new above function), this should not matter. https://github.com/lsst/daf_butler/commit/ac5032c6361cfb3a6f8d1dd78cf2b52fb6baa3e3
• Gen3 has FitsExposureFormatter, which might be able to call that new function, if we put it in a good location.
• A later addition would be to add an atmospheric refraction correction to the above model, so that the matcher does not have to account for it.
• We can directly compute this, since we know the position and sky conditions of the exposure.
• This work should be a separate ticket, once the above code is stable.

### How to test this?

We can test by looking at the initial matcher output for all sensors on a single visit. The matcher reports how far it had to move at each iteration; if we scatter plot that for all sensors in one visit, they should cluster pretty tightly if we've gotten the input CameraGeom correct. If the distortion model is not good, there will be much larger scatter and/or more outliers.

Jim suggests writing regression tests using testdata_* for the various obs packages, testing pixelToSky and skyToPixel for e.g. the ccd bounding boxes.

## ip_isr changes

ISR has an optional addDistortionModel which appears to only be turned on for HSC. This uses afw's makeDistortedTanWcs(). Although Russell also wrote this, I think we'll be much better off just using makeSkyWcs as above: it should be much simpler than what goes on in makeDistortedTanWcs.

• Disable doAddDistortionModel for all cameras, so that our above work does not get overwritten during ISR.
• Deprecate and remove the relevant code from ISR task: there is no reason to do this in ISR instead of doing it as part of the reading of the raw exposure data.

## AstrometryTask changes

Per discussion on slack[1] we will start by not writing a new fitter, but just testing how close the above initial WCS fixes get us to what we want. From there, we can decide how far we want to proceed with a new fitter.

• New fitter: FitWithFixedCameraGeomTask?
• This new fitter might not have to iterate at all. If it is just fitting crval and cdMatrix, that might be a single step: match to refcat, compute crval (still at boresight), compute cdMatrix, done.
• We want to keep crval at the boresight so that crpix is still defined at the boresight, so that the cameraGeom is unchanged. We will have to correct crval to account for any offset relative to the matched sources.
• The fitter can use a pure linear regression, fitting cdMatrix + crval offset.
• We then set AstrometryConfig.wcsFitter to this and everything else should be the same?
• While we're doing this, we might want to tidy up the the except Exception as e line in astrometry.py to not catch everything?
• If we discover problems with getting enough refcat matches (e.g. due to too large of boresight offsets), we should just increase the padding for loadPixelBox.

### How to test this?

Lauren is identifying test data for this process, which will be used by e.g. ci_hsc or test_processCcd. We can write unit tests of the fitting code itself with some existing testdata sets, and/or trivial fake points+fake refcats, just to check that the fitter is behaving correctly.

Show
John Parejko added a comment - Design/change sketch for AstrometryTask updates I think that the below encompasses all of the necessary work to implement the new single frame astrometry design. John Swinbank : is this sufficient to close out this ticket? I'm filing a new ticket to start the work on the first part. Reading raw Exposure changes This work is going to go into the Gen2 butler world. Converting it to gen3 butler will be a separate job, but could be made easier by keeping all of the relevant changes in free functions that we can just use directly. The initial SkyWcs is assigned in exposureFromImage (called from _standardizeExposure, called from std_raw, so this all happens when the raw is loaded by the butler via std_raw). exposureFromImage first calls makeSkyWcs(metadata), and only later builds the VisitInfo. We will have to reorder that so that VisitInfo gets built first, then makeskyWcs(pixelsToFieldAngle, orientation, flipX, boresight) is called later, using the camerageGeom to specify pixelsToFieldAngle and flipX, and the VisitInfo to specify orientation and boresight. Need to write a separate function to call from exposureFromImage (createInitialSkyWcs(visitInfo, detector)?) that looks similar to lsstCamMapper.getWcsFromDetector. We may have a problem due to when metadata gets stripped. Jim mentioned this gen3 commit as an example. However, if we don't actually need the FITS WCS metadata at all (which we won't with the new above function), this should not matter. https://github.com/lsst/daf_butler/commit/ac5032c6361cfb3a6f8d1dd78cf2b52fb6baa3e3 Gen3 has FitsExposureFormatter, which might be able to call that new function, if we put it in a good location. A later addition would be to add an atmospheric refraction correction to the above model, so that the matcher does not have to account for it. We can directly compute this, since we know the position and sky conditions of the exposure. This work should be a separate ticket, once the above code is stable. How to test this? We can test by looking at the initial matcher output for all sensors on a single visit. The matcher reports how far it had to move at each iteration; if we scatter plot that for all sensors in one visit, they should cluster pretty tightly if we've gotten the input CameraGeom correct. If the distortion model is not good, there will be much larger scatter and/or more outliers. Jim suggests writing regression tests using testdata_* for the various obs packages, testing pixelToSky and skyToPixel for e.g. the ccd bounding boxes. ip_isr changes ISR has an optional addDistortionModel which appears to only be turned on for HSC. This uses afw's makeDistortedTanWcs(). Although Russell also wrote this, I think we'll be much better off just using makeSkyWcs as above: it should be much simpler than what goes on in makeDistortedTanWcs. Disable doAddDistortionModel for all cameras, so that our above work does not get overwritten during ISR. Deprecate and remove the relevant code from ISR task: there is no reason to do this in ISR instead of doing it as part of the reading of the raw exposure data. AstrometryTask changes Per discussion on slack [1] we will start by not writing a new fitter, but just testing how close the above initial WCS fixes get us to what we want. From there, we can decide how far we want to proceed with a new fitter. 1: https://lsstc.slack.com/archives/GK4LKTY9W/p15604315110283000 New fitter: FitWithFixedCameraGeomTask? This new fitter might not have to iterate at all. If it is just fitting crval and cdMatrix, that might be a single step: match to refcat, compute crval (still at boresight), compute cdMatrix, done. We want to keep crval at the boresight so that crpix is still defined at the boresight, so that the cameraGeom is unchanged. We will have to correct crval to account for any offset relative to the matched sources. The fitter can use a pure linear regression, fitting cdMatrix + crval offset. We then set AstrometryConfig.wcsFitter to this and everything else should be the same? While we're doing this, we might want to tidy up the the except Exception as e line in astrometry.py to not catch everything? If we discover problems with getting enough refcat matches (e.g. due to too large of boresight offsets), we should just increase the padding for loadPixelBox. How to test this? Lauren is identifying test data for this process, which will be used by e.g. ci_hsc or test_processCcd. We can write unit tests of the fitting code itself with some existing testdata sets, and/or trivial fake points+fake refcats, just to check that the fitter is behaving correctly.
Hide
John Parejko added a comment -

Jim Bosch approves what was written above.

Show
John Parejko added a comment - Jim Bosch approves what was written above.

#### People

Assignee:
Chris Morrison
Reporter:
John Swinbank
Watchers:
Chris Morrison, Jim Bosch, John Parejko, John Swinbank, Lauren MacArthur, Tim Jenness
Votes:
0 Vote for this issue
Watchers:
6 Start watching this issue

#### Dates

Created:
Updated:
Resolved: