Manipulating WCS and CelestialCoordinates to place a location on a WFIRST SCA

[josePhoenix]

Having experimented with demo13.py and knowing a bit about the WFIRST coordinate systems, I thought I would try to simulate a piece of real sky that I have a catalog for. However, I'm not sure how to place the location I want on a particular SCA.

Supplying the coordinate I want as the world_pos= argument to wfirst.getWCS() places the target in the center of the focal plane array (so, not on any of the SCAs). I tried to work out the correct direction and angle to offset the coordinate using the stuff in galsim.wfirst.wfirst_wcs to no avail.

~I'm also having trouble understanding the plot I get from this code to show the SCA centers as projected on the sky. It looks like there may be a parity flip somewhere (comparing to this plot: https://github.com/GalSim-developers/GalSim/wiki/GalSim-WFIRST-module-diagrams).~

Never mind, RA should increase to the left...

from matplotlib import pyplot as plt
import galsim
from galsim import wfirst
from galsim.wfirst import wfirst_wcs

import datetime
vernal_equinox_2025 = datetime.datetime(2025,3,20,9,2,0)
subframe_ra, subframe_dec = (11.721605694716114, 85.277622661476954)
targ_pos = galsim.CelestialCoord(
    ra=subframe_ra * galsim.degrees,
    dec=subframe_dec * galsim.degrees,
)
wcs_lookup = wfirst.getWCS(world_pos=targ_pos, date=vernal_equinox_2025)
sca_ra_centers, sca_dec_centers = [], []
for i in range(1, 19):
    center = wcs_lookup[i].toWorld(galsim.PositionD(wfirst.n_pix/2, wfirst.n_pix/2))
    ra, dec = center.ra / galsim.degrees, center.dec / galsim.degrees
    plt.scatter(ra, dec)
    plt.text(ra, dec, str(i))
plt.scatter(subframe_ra, subframe_dec)

[rmjarvis]

I suspect this is the Cycle 5 -> Cycle 6 update that is currently being implemented by @rmandelb. Cross referencing Issue #675.

[rmandelb]

Hi - @rmjarvis is right that the WCS is going to updated due to some Cycle 5 -> 6 updates. You can find a preview on the branch that he mentioned. There are a few updates:

• it seems there might be a slight overall offset of the entire focal plane compared to what there should be
• in cycle 5->6, they swapped around 2 rows of SCA labels

I am waiting to hear from some collaborators to confirm that the new code is OK, but after that, it will be updated.

However, it seems you had two somewhat separate issues. It seems like the parity flip issue might not have been real (from your update in the initial issue) but it sounds like your first concern was a difficulty in making a correspondence between some position that you want to observe, and the actual SCA it lives on? Have I understood this correctly?

As outlined in the docstring for galsim.wfirst.getWCS(), the world_pos argument does indeed correspond to the center of the focal plane as on the diagram here: https://github.com/GalSim-developers/GalSim/wiki/GalSim-WFIRST-module-diagrams (marked as a plus sign inside a circle) Currently, there is no option to specify that world_pos corresponds to something else, like the center of one of the SCAs. In your shoes, I would probably do the following:

1. Choose some central position in your catalog as the center of the focal plane array, and get the full set of WCSs for all the SCAs.
2. For each object in the catalog, use its position along with the galsim.wfirst.findSCA routine to figure out which SCA it lives on (if any; it might of course be in a gap or off the edge).
3. Then use the WCS for the SCA that was identified in step (2) to get GalSim to draw the object in the right image, with the right location and orientation.

If your catalog covers a decent area, this should be fine. If it's a tiny area, like less than a single SCA, this is probably a really clunky way to do it as you have too high a chance to end up with most of the objects in a gap.

Does that make sense? Sorry, I may not have fully grasped what you are trying to do, but I'm happy to discuss further if needed.

[josePhoenix]

I think you've understood my issue correctly.

If your catalog covers a decent area, this should be fine. If it's a tiny area, like less than a single SCA, this is probably a really clunky way to do it as you have too high a chance to end up with most of the objects in a gap.

That's exactly my problem. I fiddled around with an offset to the world_pos to place my target on SCA 1, but couldn't work out a way to calculate the appropriate offset from the geometry information in wfirst_wcs. It seems like it should be possible, at least as long as the PA is fixed, but I stared at it for a while yesterday without success.

[rmandelb]

Do you actually care about having the right RA/dec? You could do a best-guess nearby position for world_pos, find the coordinates that actually live on some SCA (pick your favorite, I guess!), then shift the catalogs RA/dec values appropriately so they live on the SCA. It seems like this should be fine if you just want to get an image to stare at. If you want to match it against some other image then that won't work.

In terms of how to better support this use case, I think the way to make it easier for a user to do this would be:

• take the chunk of code in getWCS() that finds the position on the sky for the center of each SCA given world_pos and the orientation of the observatory, and turn it into a callable routine that is exposed to the user
• getWCS() can then call that routine when needed
• users can call that routine in order to back out where SCAs would live given some world_pos and orientation

Does that sound like it would be helpful? I could add it to the feature requests in #675 and hopefully address them in the coming weeks.

[josePhoenix]

That does sound helpful. I feel like I should be able to do this with some trigonometry, but an analytic solution eludes me. This is what I've come up with for now:

def center_location_on_sca(ra, dec, sca, date=vernal_equinox_2025):
    targ_pos = galsim.CelestialCoord(
        ra=ra * galsim.degrees,
        dec=dec * galsim.degrees,
    )
    best_pa = wfirst.bestPA(targ_pos, vernal_equinox_2025)
    def _optimize(coords):
        fpa_ra, fpa_dec = coords
        fpa_pos = galsim.CelestialCoord(
            ra=fpa_ra * galsim.degrees,
            dec=fpa_dec * galsim.degrees,
        )
        wcs_lookup = wfirst.getWCS(world_pos=fpa_pos, date=date, SCAs=sca)
        sca_pos = wcs_lookup[sca].toWorld(galsim.PositionD(wfirst.n_pix/2, wfirst.n_pix/2))
        delta = targ_pos.distanceTo(sca_pos) / galsim.radians
        return delta
    result = optimize.fmin(_optimize, [ra, dec])
    return result

[rmjarvis]

I think this issue can probably be closed. We never implemented the above suggestion, but it seems like Jose figured out something that worked for what they were doing at the time. Please feel free to reopen if someone thinks this or a similar use case would benefit from some additional helper code in GalSim.