add match_source_position option for setup_sim_to_match_file

[mperrin]

The function setup_sim_to_match_file has been setting the detector_position attribute to be the center of the detector (or center of subarray, etc, depending on the aperture name). The target may not necessarily be located there, particularly for dithered images.

This PR adds a new keyword match_source_position to that function, which uses the FITS header TARG_RA, TARG_DEC and the WCS keywords to figure out the pixel coordinates where the target is expected to be, and it sets up the simulation instance to have detector_position set to that, rounded to integer pixels. (It also attempts to set sub pixel position using source_offset, even though this probably won't be precise enough in many cases)

Questions:

1. Should this maybe be turned on by default? For now I have left it off by default, for back compatibility so that existing code doesn't get unexpected changes of behaviors. But there's a case to me made that this might as well be on by default in general. ?
2. Is trying to set sub pixel position with source offsets helpful, or confusing and unnecessary? Should there be a separate option to toggle this part on and off? This can only be as good as the WCS, which is not in general always going to be precise enough at a sub pixel level.

Example:

In [3]:  miri = stpsf.setup_sim_to_match_file('jw01193062001_02101_00001_mirimage
      ⋮ _cal.fits', match_source_position=True)

Setting up sim to match jw01193062001_02101_00001_mirimage_cal.fits
iterating query, tdelta=3.0

MAST OPD query around UTC: 2022-11-13T06:57:01.514
                        MJD: 59896.28960085648

OPD immediately preceding the given datetime:
	URI:	 mast:JWST/product/R2022111308-NRCA3_FP1-1.fits
	Date (MJD):	 59896.0177
	Delta time:	 -0.2719 days

OPD immediately following the given datetime:
	URI:	 mast:JWST/product/R2022111503-NRCA3_FP1-1.fits
	Date (MJD):	 59898.0069
	Delta time:	 1.7173 days
User requested choosing OPD time closest in time to 2022-11-13T06:57:01.514, which is R2022111308-NRCA3_FP1-1.fits, delta time -0.272 days
Importing and format-converting OPD from /Users/mperrin/software/webbpsf-data/MAST_JWST_WSS_OPDs/R2022111308-NRCA3_FP1-1.fits
Backing out SI WFE and OTE field dependence at the WF sensing field point (NRCA3_FP1)
Sensing inst model using apername NRCA3_FP1
Using sensing field point SI WFE model from file wss_target_phase_fp1.fits

Configured simulation instrument for:
    Instrument: MIRI
    Filter: F2550W
    Detector: MIRIM
    Apername: MIRIM_BRIGHTSKY
    Det. Pos.: (256, 256) in subarray
    Image plane mask: None
    Pupil plane mask: None

new part starts here:

    Attempting to match target position precisely...
    Target RA, Dec from file header: 22h42m22.087s -29d21m39.77s
    Based on WCS that is at pixel coords [326.027730802757, 196.1974238785417]
    Setting Det. Pos.: (326, 196) in subarray
    Setting subpixel offset: [0.0277308  0.19742388] pixels = 0.003, 0.022 arcsec.

[marcio-melendez]

This is a bit tricky with large mosaic observations because the target RA/DEC is the same for all the tiles but some tiles may be one full detector away from the target, (X/YOFFSET from the dither information header). So, when trying a cal_fit from a large mosaic I get the following error:

ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 psf = stpsf.setup_sim_to_match_file('./jw02739010001_02103_00001_nrcalong_cal.fits', match_source_position=True)

File ~/stpsf_setup_sim_source_pos/stpsf/match_data.py:142, in setup_sim_to_match_file(filename_or_HDUList, verbose, plot, choice, match_source_position)
    139 targ_pos_integer_part = np.asarray(np.round(targ_pos), int)
    140 targ_pos_subpix_part = targ_pos - targ_pos_integer_part
--> 142 inst.detector_position = targ_pos_integer_part
    143 inst.options['source_offset_x'] = targ_pos_subpix_part[0] * inst.pixelscale
    144 inst.options['source_offset_y'] = targ_pos_subpix_part[1] * inst.pixelscale

File ~/stpsf_setup_sim_source_pos/stpsf/stpsf_core.py:318, in SpaceTelescopeInstrument.detector_position(self, position)
    316     raise ValueError('Detector pixel coordinates must be pairs of nonnegative numbers, not {}'.format(position))
    317 if x < 0 or y < 0:
--> 318     raise ValueError('Detector pixel coordinates must be nonnegative integers')
    319 if isinstance(self._detector_npixels, tuple):
    320     # A tuple has been provided for a non-square detector with different Y and X dimensions
    321     det_npix_y, det_npix_x = self._detector_npixels

ValueError: Detector pixel coordinates must be nonnegative integers

In addition, when trying to use this new functionality with the wavefront sensing level 3 products, I get the following error about WCS coordinates not found.

File ~/anaconda3/envs/py311_2025/lib/python3.11/site-packages/stdatamodels/properties.py:349, in ObjectNode.__getattr__(self, attr)
    348 try:
--> 349     val = self._instance[attr]
    350 except KeyError as err:

KeyError: 'wcs'

The above exception was the direct cause of the following exception:

AttributeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 psf = stpsf.setup_sim_to_match_file('./jw04503-o091_t019_nircam_f212n-wlm8-nrca1_wfscmb-04.fits', match_source_position=True)

File ~/stpsf_setup_sim_source_pos/stpsf/match_data.py:138, in setup_sim_to_match_file(filename_or_HDUList, verbose, plot, choice, match_source_position)
    134 targ_coords = astropy.coordinates.SkyCoord(model.meta.target.ra,
    135                             model.meta.target.dec,
    136                             frame='icrs', unit=u.deg)
    137 # where does that show up in the FITS file
--> 138 targ_pos = model.meta.wcs.world_to_pixel(targ_coords)
    139 targ_pos_integer_part = np.asarray(np.round(targ_pos), int)
    140 targ_pos_subpix_part = targ_pos - targ_pos_integer_part

File ~/anaconda3/envs/py311_2025/lib/python3.11/site-packages/stdatamodels/properties.py:352, in ObjectNode.__getattr__(self, attr)
    350 except KeyError as err:
    351     if schema == {}:
--> 352         raise AttributeError(f"No attribute '{attr}'") from err
    354     val = _make_default(attr, schema, self._ctx)
    355     if val is not None:

AttributeError: No attribute 'wcs'

[mperrin]

Good catches @marcio-melendez, thanks. I will work at adding some extra logic to handle these cases...

Do you think the preferable behavior in such cases would be:

1. To raise a ValueError or RuntimeError that stops the code, with an informative message to the user to let them know there is an issue, or
2. To catch the exception, and just skip trying to set the position, and leave the position at the current default center coordinates?

We can easily make it do either, so we just need to pick which of those would be preferable to users.

[marcio-melendez]

Good catches @marcio-melendez, thanks. I will work at adding some extra logic to handle these cases...

Do you think the preferable behavior in such cases would be:

1. To raise a ValueError or RuntimeError that stops the code, with an informative message to the user to let them know there is an issue, or
2. To catch the exception, and just skip trying to set the position, and leave the position at the current default center coordinates?

We can easily make it do either, so we just need to pick which of those would be preferable to users.

I think option 2. Also, we'll need to update the documentation page for the "setup simulation from observation" to capture this new option, we could add that change to this PR.

[marcio-melendez]

I am getting an name error for this variable, asdasd

[mperrin]

ha that's clearly some typo junk... my mistake! I don't remember what went on with that commit...

[marcio-melendez]

I removed the extra word that I mentioned before and I was about to add that typo change to this commit but then I noticed this "issue".
In cases with invalid coordinates there is a warning when trying to calculate the integer part of the position. Later there is a raise ValueError because of that. Maybe there is no need to try to calculate the new position if the coordinates are invalid, just stop there and raise the Value Error from that point. Otherwise we'll have unnecessary warning messages.
RuntimeWarning: invalid value encountered in cast
targ_pos_integer_part = np.asarray(np.round(targ_pos), int)