Uploaded image for project: 'Data Management'
  1. Data Management
  2. DM-14932

Add utility functions for creating SkyWcss from boresight/rotator + cameraGeom

    XMLWordPrintable

Details

    • Story
    • Status: Done
    • Resolution: Done
    • None
    • afw
    • 2
    • AP F18-2
    • Alert Production

    Description

      Given an accurate boresight and rotator angle from a camera, we should be able to construct our best initial estimate of the WCS for a sensor by combining that information with our own camera geometry, as this will probably produce better results that relying on any WCS provided directly by the camera.  In fact, for LSST, it's not clear that there even will be a sensor WCS provided directly by the camera, because LSST raw data has amps in different HDUs and no full-sensor HDU that could unproblematically hold a full-sensor WCS.

      To do without using astshim code directly, we should provide utility functions in afw with something like the following signatures:

      TransformPoint2ToSpherePoint makeIntermediateWorldCoordsToSky(
          SpherePoint const & position,  // (ra, dec) that corresponds to FIELD_ANGLE origin
          Angle rotation // rotation about FIELD_ANGLE origin
      );
       
      std::shared_ptr<SkyWcs> makeSkyWcs(
          TransformPoint2ToPoint2 const & pixToIwc,
          TransformPoint2ToSpherePoint const & iwcToSky
      );

      Attachments

        Activity

           I believe that this function works (it might be easier once e.g. makeIntermediateWorldCoordsToSky is available).

          The API is higher than the examples given above, but I think that it needs to be supported as part of this ticket (the details of the API could change, of course, but I deliberately chose to pass a Point2D and Angle and would be a little unhappy if this were changed as that seems to me to be the natural inputs for the user).

          When DM-14980 is resolved the API will not need to include the camera.

           

          import astshim
          import lsst.afw.cameraGeom as cameraGeom
          import numpy as np                      # only for pi!
           
          def getWcsFromDetector(camera, detector, boresight, rotation=0*afwGeom.degrees):
              """Given a camera, detector and (boresight, rotation), return that detector's WCS
           
                  Parameters
                  ----------
                  camera: `lsst.afw.cameraGeom.Camera`  The camera containing the detector
                  detector: `lsst.afw.cameraGeom.Detector`  A detector in a camera
                  boresight: `lsst.afw.geom.Point2D`  The boresight of the observation
                  rotation: `lsst.afw.geom.Angle` The rotation angle of the camera
           
              The rotation is "rotskypos", the angle of sky relative to camera coordinates
              (from North over East)
              """
              trans = camera.getTransform(detector.makeCameraSys(cameraGeom.PIXELS),
                                          detector.makeCameraSys(cameraGeom.FIELD_ANGLE))    
              polyMap = trans.getMapping()
              radToDeg = astshim.ZoomMap(2, 180/np.pi) # convert from radians to degrees
              polyMap = polyMap.then(radToDeg)
           
              pixelFrame = astshim.Frame(2, "Domain=PIXELS")
              iwcFrame = astshim.Frame(2, "Domain=IWC")
           
              frameDict = astshim.FrameDict(pixelFrame, polyMap, iwcFrame)
           
              crpix = afwGeom.PointD(0, 0)
              crval = afwGeom.SpherePoint(*boresight, afwGeom.degrees)
              cd = afwGeom.makeCdMatrix(1.0*afwGeom.degrees, -(rotation + 90*afwGeom.degrees), True)
              iwcToSkyWcs = afwGeom.makeSkyWcs(crpix, crval, cd)
           
              iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
              skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
           
              frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
           
              wcs = afwGeom.SkyWcs(frameDict)
              
              return wcs
          

           

          rhl Robert Lupton added a comment -  I believe that this function works (it might be easier once e.g. makeIntermediateWorldCoordsToSky is available). The API is higher than the examples given above, but I think that it needs to be supported as part of this ticket (the details of the API could change, of course, but I deliberately chose to pass a Point2D and Angle and would be a little unhappy if this were changed as that seems to me to be the natural inputs for the user ). When DM-14980 is resolved the API will not need to include the camera.   import astshim import lsst.afw.cameraGeom as cameraGeom import numpy as np # only for pi!   def getWcsFromDetector(camera, detector, boresight, rotation = 0 * afwGeom.degrees): """Given a camera, detector and (boresight, rotation), return that detector's WCS   Parameters - - - - - - - - - - camera: `lsst.afw.cameraGeom.Camera` The camera containing the detector detector: `lsst.afw.cameraGeom.Detector` A detector in a camera boresight: `lsst.afw.geom.Point2D` The boresight of the observation rotation: `lsst.afw.geom.Angle` The rotation angle of the camera   The rotation is "rotskypos" , the angle of sky relative to camera coordinates ( from North over East) """ trans = camera.getTransform(detector.makeCameraSys(cameraGeom.PIXELS), detector.makeCameraSys(cameraGeom.FIELD_ANGLE)) polyMap = trans.getMapping() radToDeg = astshim.ZoomMap( 2 , 180 / np.pi) # convert from radians to degrees polyMap = polyMap.then(radToDeg)   pixelFrame = astshim.Frame( 2 , "Domain=PIXELS" ) iwcFrame = astshim.Frame( 2 , "Domain=IWC" )   frameDict = astshim.FrameDict(pixelFrame, polyMap, iwcFrame)   crpix = afwGeom.PointD( 0 , 0 ) crval = afwGeom.SpherePoint( * boresight, afwGeom.degrees) cd = afwGeom.makeCdMatrix( 1.0 * afwGeom.degrees, - (rotation + 90 * afwGeom.degrees), True ) iwcToSkyWcs = afwGeom.makeSkyWcs(crpix, crval, cd)   iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping( "PIXELS" , "SKY" ) skyFrame = iwcToSkyWcs.getFrameDict().getFrame( "SKY" )   frameDict.addFrame( "IWC" , iwcToSkyMap, skyFrame)   wcs = afwGeom.SkyWcs(frameDict) return wcs  
          rowen Russell Owen added a comment - - edited

          I added one new function: makeSkyWcs(pixelsToFieldAngle, orientation, flipX, crval, projection="TAN") rather than the two requested functions. I did this because it is not enough to have the IWC to sky transform, one must also construct an appropriate ast::SkyFrame with the SkyRef and SkyRefIs properties set, and I wanted to avoid users having to deal with that or with FrameDicts.

          Note that "orientation" is with respect to the focal plane, not CCD x, y. This can be helpful for cameras such as HSC that have a mix of CCD rotations, because the same angle can be used for all CCDs. However, it does require some care on the part of the user.

          rowen Russell Owen added a comment - - edited I added one new function: makeSkyWcs(pixelsToFieldAngle, orientation, flipX, crval, projection="TAN") rather than the two requested functions. I did this because it is not enough to have the IWC to sky transform, one must also construct an appropriate ast::SkyFrame with the SkyRef and SkyRefIs properties set, and I wanted to avoid users having to deal with that or with FrameDicts . Note that "orientation" is with respect to the focal plane, not CCD x, y. This can be helpful for cameras such as HSC that have a mix of CCD rotations, because the same angle can be used for all CCDs. However, it does require some care on the part of the user.
          rowen Russell Owen added a comment -

          For the record: the FITS intermediate world coordinates frame (IWC) is the same as field angle, except as follows:

          • IWC is in degrees, whereas field angle is in radians
          • flipX affects IWC, but not field angle
          • orientation affects IWC, but not field angle
          rowen Russell Owen added a comment - For the record: the FITS intermediate world coordinates frame (IWC) is the same as field angle, except as follows: IWC is in degrees, whereas field angle is in radians flipX affects IWC, but not field angle orientation affects IWC, but not field angle
          jbosch Jim Bosch added a comment -

          Looks good.  Some comments about tests and names on the PR.

          I'll also note that I do think that the use of SpherePoint for the boresight position here is preferable to the Point2d in rhl's proposed signature above; we should use Angle for angles, and SpherePoint is the way to do that for coordinate tuples.

          But the functionality to go straight from a Detector to a SkyWcs in rhl's proposed signature definitely is desirable, even if it's just a trivial forward to the function implemented here that operates on the Transform.  Given that that signature would be changing soon anyway on DM-14980, I'm happy to just add it when I take that ticket on, but I also wouldn't be opposed to rowen adding it here.  I'll also admit that I haven't thought yet about whether it'd be problematic to have that function here from a circular dependency standpoint (I don't know which way the afw.geom<->afw.cameraGeom dependency runs now that we've moved some stuff into lsst.geom).  Regardless, we'll want that signature somewhere.

          jbosch Jim Bosch added a comment - Looks good.  Some comments about tests and names on the PR. I'll also note that I do think that the use of SpherePoint for the boresight position here is preferable to the Point2d in  rhl 's proposed signature above; we should use Angle for angles, and SpherePoint is the way to do that for coordinate tuples. But the functionality to go straight from a Detector to a SkyWcs in rhl 's proposed signature definitely is desirable, even if it's just a trivial forward to the function implemented here that operates on the Transform .  Given that that signature would be changing soon anyway on DM-14980 , I'm happy to just add it when I take that ticket on, but I also wouldn't be opposed to rowen adding it here.  I'll also admit that I haven't thought yet about whether it'd be problematic to have that function here from a circular dependency standpoint (I don't know which way the afw.geom<->afw.cameraGeom dependency runs now that we've moved some stuff into lsst.geom).  Regardless, we'll want that signature somewhere.
          rowen Russell Owen added a comment -

          Thank you for the helpful review jbosch

          rowen Russell Owen added a comment - Thank you for the helpful review jbosch

          People

            rowen Russell Owen
            jbosch Jim Bosch
            Jim Bosch
            Jim Bosch, John Parejko, Robert Lupton, Russell Owen
            Votes:
            0 Vote for this issue
            Watchers:
            4 Start watching this issue

            Dates

              Created:
              Updated:
              Resolved:

              Jenkins

                No builds found.