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

Error in computation of MIRI LRS fixed slit S_REGION footprint #9072

Open
stscijgbot-jp opened this issue Jan 13, 2025 · 49 comments · May be fixed by #9193
Open

Error in computation of MIRI LRS fixed slit S_REGION footprint #9072

stscijgbot-jp opened this issue Jan 13, 2025 · 49 comments · May be fixed by #9193

Comments

@stscijgbot-jp
Copy link
Collaborator

Issue JP-3848 was created on JIRA by David Law:

There may be an existing ticket for this, but if so I can't find it.

LRS SLIT data is giving incorrect S_REGION results that don't have the right size for the real LRS slit.  After some poking at the issue:

The initial S_REGION in the uncal.fits files looks reasonable, the issue is with the updated S_REGION computed during the assign_wcs step.

The root issue seems to be in the gwcs.footprint code.  This is passed in the LRS slit bounding box, which covers the entire ~ 400 row high region of the detector onto which LRS can disperse light.  This is then mapped back to the sky using the LRS distortion solution.  However, since the LRS distortion solution is tied into the MIRI imaging distortion solution at a fixed wavelength, the bounding box actually maps to a line in the sky in RA/DEC coordinates that's tracing the middle of the slit.  This line is then projected to a box based on the max/min RA/DEC coordinates, and is thus perfectly aligned with RA/DEC with a size in each direction that depends on the spacecraft roll angle.

The S_REGION reported after the assign_wcs step is thus variable size based on the roll angle, and has an incorrect rotation.

For the LRS slit, the bounding box will not be useful to describe the slit footprint and some other method will presumably need to be used instead (e.g., the SIAF aperture).

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Assigning to Jane Morrison since I know you're already working on it (and potentially solved it already on some other ticket?).

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Example with a nearly square S_REGIONS slit footprint (since roll angle is near 45 degrees): jw03168004001_04103_00001_mirimage

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

For reference, s_region.png attached shows the computed S_REGION footprint translated back into detector pixel coordinates for two different roll angle cases.  This should always look like the illuminated region of the slit, but can actually have dramatically different shapes.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Also tagging Nadia Dencheva FYI.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law or Nadia Dencheva 

Where is the S_REGION calculated in a uncal file ? 
By that I mean what routine is used to determine this value in an uncal file ?

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Nadia Dencheva on JIRA:

Jane Morrison The routine for the uncal file uses the SIAF and does not apply to cal or s2d files. It's in set_telescope_pointing.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

Below is a slack conversion. I have pulled it in this ticket to track the history of discussions. 
dlawdlaw:
A few quick thoughts: * First, I don't think it will be possible to use the WCS reference files in CRDS to derive the LRS slit footprint.  The bounding box is in detector pixels and just describes the region of the detector that can have light dispersed over it, and there's nothing about the WCS than intrinsically knows where the slit corners are.  Since there's spatial/spectral degeneracy in the slit, the WCS maps everything back to a single line at constant wavelength and doesn't know anything about the slit width.

  • In contrast, the MIRIM_SLIT aperture in SIAF is of AperType=SLIT meaning that it's defined in V2V3 space (or strictly in Ideal coordinates that map directly to V2V3 space with a simple rotation).  I.e., there's no distortion to worry about- the values are already mapped through detector distortion to V2V3.  As such, all that needs to be done to map the SIAF entry to RA/DEC is to apply the telescope attitude matrix (and presumably the velocity aberration correction).
  • I don't think that we can query the SIAF directly in assign_wcs because pysiaf is not a dependency of jwst.  It's used in SDP, but not cal pipeline.

So, what to do?  In this case, it may be as simple as NOT recomputing S_REGION in assign_wcs and just keeping the one populated by set_telescope_pointing in SDP.  Looking at those values, they match what I'd expect for the true slit footprint on the sky (based on passing the known V2V3 corners into the v2v3 to world transform in the data model after assign_wcs has been run) to within 0.1 milliarcsec.  I believe that's consistent with a lack of differential velocity aberration correction in the S_REGION populated by SDP, and should be irrelevant?
_______________________
denchevadencheva  9:11 AM

Aren't the corners of the slit in the reference file? Can these be used with the WCS to compute a footprint?
 _______________________
dlawdlaw  9:11 AM

I checked and no, the corners of the slit are not in the reference file.  We could certainly add them though.
(and if we added the corners we could certainly use those with the WCS to compute a footprint)
___________________________

 
denchevadencheva  9:23 AM

Then the best solution would be adding them to the reference file and for the time being just not recomputing S_REGIOON in assign_wcs and resample.
!https://ca.slack-edge.com/EQN4UL0NM-WRX8GJTEZ-g0c69db50a80-48!
morrisonmorrison  9:26 AM

What is the difference to what values would be added to the reference file and what is in the SIAF now. Would we be off if we just did not recalculate the S_REGION ?
!https://ca.slack-edge.com/EQN4UL0NM-WRZE6U0CE-99f960a1f754-48!
denchevadencheva  9:27 AM

We'll be off but not by much since the slit is small and the distortion over that is not that big
!https://ca.slack-edge.com/EQN4UL0NM-WRKPSAN1H-500ccb3b57e8-48!
dlawdlaw  9:27 AM

Strictly speaking we're off if we don't recalculate because the residual velocity aberration isn't taken into account, even if it is absolutely tiny.
!https://ca.slack-edge.com/EQN4UL0NM-WRX8GJTEZ-g0c69db50a80-48!
morrisonmorrison  9:28 AM

@dlaw who would be the person to contact about getting the values of the slit corners in the reference file. And would we then want to apply a residual velocity aberration correction to those values ?
----New
!https://ca.slack-edge.com/EQN4UL0NM-WRKPSAN1H-500ccb3b57e8-48!
dlawdlaw  9:30 AM

@morrison I'd add that information in the LRS reference file myself, in v2v3 coordinates.  All we'd need to do then is apply the distortion model to it: model.meta.wcs.transfom('v2v3', 'world', v2, v3, lam)  where the lambda value is probably necessary for format but irrelevant.It's probably worth taking a quick look at some other slit-type modes though just to make sure that those don't need attention too.
9:31
dlaw
(since the va correction is part of the distortion model chain internally)

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law for testing if you can send me the values you intend to put in the reference file then I can use those in code to make sure I am using them correctly and getting reasonable S_REGION values.  Or if you have the updated reference file rather soon - send me the I can test with that (before it is in CRDS).

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

v2,v3= array([-412.54598854, -412.50298855, -417.21103301, -417.25301145]), array([-401.01798362, -400.50498252, -400.11702824, -400.63101748])

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I took a quick look at a NIRSpec slit case and there doesn't seem to be the same kind of rotation issue, but there is an unexpected mismatch between the S_REGION information in the cal.fits files and what I'd expect based on the pysiaf slit aperture.

See nirspec.png attached for NRS_S200A1_SLIT in jw02072004001_03102_00001_nrs1_rate

The slit as described by S_REGION is somewhat longer than defined by the pysiaf aperture.  Is this a known thing Christian Hayes ?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Christian Hayes on JIRA:

David Law yes this is known for NIRSpec.  The cal SREGION should be 10% larger than the pysiaf aperture.  When performing the 2D extraction of NIRSpec fixed slits there is a 5%, buffer to either side of the extracted slit region in the cross-dispersion direction to avoid having the edge of the slit right at the edge of the 2D cal spectra (and possibly losing a pixel or two at the edge of the slit).  This means that the equivalent on sky area of the cal spectra is slightly larger than the actual slit size.

I'm not sure if the S_REGION is used for anything (other than potentially by users).  If it isn't this could be updated to only give the illuminated slit region, through it would require keeping track of that buffer size (because it differs between MOS and FS and I believe they use the same function for producing the S_REGION).

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law what about the slitless case ?

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

LRS slitless S_REGIONS look fine to me; in that case the region really is defined by the pixel bounding box values, so the current computation works fine.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

First tests look reasonable. 
I compared the s_region that is in the rate file (filled in by set_telescope_pointing) with the s_region from the current method and the differences are

Ra corner differences in arc seconds

3.3988032000031865, 3.8090916000001585, -3.369920399998705, -3.78011880000102
Dec corner differences in arc seconds 

-0.0546120000038286, -3.4498475999996003, 0.027687599988723832, 3.421508399995332

Using the corners of the V2, V3 of the slit that David gave me and using the model transform of v2,v3 to world (which accounts for the velocity aberration correction) I get differences with the original corners determined from set_telescope_pointing:

Ra difference in arc seconds: 
9.91826460960965e-05, 0.00012786427916466891, -0.00010101971099629736, -0.0001275188338922817

Dec difference in arc seconds:

-0.00012842505299204277, -0.00010203120268670318, 0.00012808504976646873, 0.0001033574761777345

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Attached a draft reference file jwst_miri_specwcs_2025jan17.fits
This has new keywords V2VERT1-4 and V3VERT1-4 in the primary header giving the V2V3 coordinates of the slit aperture

A few caveats:

  • These new keywords aren't reflected in the SpecwcsModel data model, but that data model already doesn't include existing important LRS header keywords such as IMX and IMY.  Likewise, the slitless specwcs reference file wouldn't have these keywords as they only make sense for the slit.  We could consider making special datamodels for the different LRS specwcs files, but given that we may already need to open this box of worms for the implementation of LRS WFSS later this year I'm inclined to punt on this for now.
  • The S_REGION update routine will have to be able to access the specwcs reference file to get these values, which will only be relevant for LRS_SLIT.
  • MIRI team will have to remember to keep these keywords updated in new deliveries of the specwcs reference file.  E.g., if the MIRI distortion model changes (expected sometime later this year) the v2v3 footprint as described in the specwcs reference file will also have to change and stay in sync with pysiaf.  (Although that's no different than the need for distortion reference files and pysiaf to stay in sync).

Tagging Sarah Kendrew and Greg Sloan so that they're aware of the above.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Nadia Dencheva on JIRA:

The MIRI LRS specwcs file is the only reference file without data model. We should write one. Currently the pipeline opens the FITS file directly. We could consider changing the format to asdf if the INS team is open to this and we can add the new keywords at that time. Using a model to generate the reference file would make it easier to remember to add the keywords.

Sarah Kendrew This should solve the issue in  INC203853.

 

@jemorrison
Copy link
Collaborator

I can create a datamodel for the MIRI LRS specwcs. @skendrew @drlaw1558 I can take a stab at converting the fits one to asdf. I would start with jwst_miri_specwcs_2025jan17.fits unless someone on the MIRI team wants to make this conversion.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law I have a draft PR that uses hardcoded V2,V3 coordinates to determine the s_region. But I have a basic question - do you need to change the lrs pipeline defined in assign_wcs util.py ? I should push up this draft pr (with the understanding that I will change the hard coded values to reading in from the reference file. ) I by pass how the S-region is updated now and it starts with the new V2,V3 slit corners. If old method was causing problems with how s_region was determine I am starting to think we need to update the lrs pipeline - just not sure if that is true - if it is not please explain why it does not work for S-region but does for mapping a pixel to sky.

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Jane Morrison The LRS team should weigh in on whether or not to convert the specwcs file to ASDF; my feeling is that it will be easiest to keep them in FITS format (Sarah Kendrew Greg Sloan ) since there may be team-based analysis tools that assume the FITS format.

We should also keep in mind that there are two different LRS specwcs reference files, one for slit and one for slitless.  We can either make one datamodel for each, or have them share a datamodel, in which case we'd need to add V2VERT/V3VERT header information for the slitless file as well.  I don't see any intrinsic problem with doing that.

We shouldn't need to make any updates to the lrs pipeline; the reason the old method didn't work is because mapping from pixel space to v2v3 first dealt with the dispersed trace wavelength information to go to the slit center line, and then through the imaging WCS to go to v2v3.  (And even if we could use the imaging WCS directly, the LRS detector bounding box is very different from the slit footprint).  In contrast, mapping from v2v3 to ra/dec just uses the later stages of the transform chain based on header attitude information, so it no longer cares about any MIRI-specific issue.

@jemorrison
Copy link
Collaborator

@drlaw1558 Ok I still have a few clarification questions.
I am still murky on why it it ok to map a pixel from detector to abl and then abl to v2v3 if it is not correct to do this for determining the s_region.

The lrs pipeline defined in assign_wcs.py miri.py is

pipeline = [(detector, dettoabl),
            (miri_focal, abltov2v3l),
            (v2v3, va_corr),
            (v2v3vacorr, teltosky),
            (world, None)]

the dettoabl is defined by lrs_xytoabl which uses a transform through the center of slit. So we still mapping the x,y through the center of slit to get to abl. Is the answer that we can do this for a single pixel because wavelength issue is not important when trying to find ra, dec. But for the s_region when we use these transforms it is incorrect but I am still struggling with how the roll angle of the telescope comes into play.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Greg Sloan on JIRA:

It would make sense to retain the current FITS table format.

@jemorrison
Copy link
Collaborator

@gcsloan Ok so we will retain the FIT table format but I will go ahead and make a datamodel for it. @drlaw1558 I think it might make sense to have a datamodel fixed slit and one for slitless since the V2,V3 corners are not needed for slitless case. It might be cleaner if we keep that separate - but if the MIRI team prefers to have the same format for both fix and slitless we can make 1 model

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Greg Sloan on JIRA:

Would it help guide us if I posted the code I use to generate the spec_wcs file? I have to set a small number of FITS keywords which should probably be part of the model.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

There's a couple reasons:

First, the current S_REGION is determined from the WCS bounding box in pixel coordinates.  This is the area over the detector for which the WCS should be defined, which is (very roughly) 50 pixels wide by 400 pixels high.  If we were to transform that to v2v3 using the imager distortion model it would be very wrong, as the actual slit opening that sees the sky is much smaller.

If we determined the S_REGION from the actual corners of the slit in pixel coordinates the imager distortion model would give the right solution, but the LRS distortion model would not since the detector -> alphabetalambda and alphabetalambda -> v2v3lambda transforms first apply the specwcs trace information to map back to the slit center coordinates, and then the imager distortion solution on those slit center coordinates.  That gives a straight line in v2v3 and in RA/DEC.  The existing S_REGION computation code assumes a rectangular region in RA/DEC (with 4 corresponding corner coordinates) rather than a line (with just two corner coordinates) as input.  If there are just two corner coordinates it assumes that must be because the other two are degenerate, corresponding to a box aligned with RA/DEC where the diagonal is the slit center (and hence the box size depends on the roll).

With the new approach we cut out all of that, and just pass in the known slit V2V3 corners and transform them to RA/DEC corners using the wavelength-independent spacecraft attitude information in the v2v3-world transform.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Jane Morrison Greg Sloan I can get behind two data models; there are some other keywords as Greg mentioned that might make sense to split between the two.  For instance, probably no need to record the slitless trace center point in the slit header information, and vice versa (right now I think the same info is copied into both files).  That also means that when LRS WFSS is introduced we'd have a third datamodel of TBD format.

@jemorrison
Copy link
Collaborator

@gcsloan I can use the existing spec wcs to set up the data models. If you plan on adding or changing the spec wcs as part of this effort then I can add/adapt the datamodel. If you do plan on changing the slitless spec wcs it would be useful to know what those changes will be and maybe if you have the time create a new fits file - then I can make sure the datamodel is what you want.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Greg Sloan on JIRA:

The only changes I know to expect are that we will soon be (1) tilting the wavelength elements slightly, and (2) varying that tilt with row position. Neither of these changes should require changes to the data model.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law Nadia Dencheva !s_region_3168.png!
FYI:
I went down a bit of a rabbit hole. I have looked at many LRS Fixed slit s_regions and they have always seemed reasonable. As I was testing this ticket I noticed some very strange results when using the develop version of the JWST pipeline from main. So I went and looked at the bulk of the MIRI LRS Fixed slit data that I have worked with and it was processed using jwst 1.16.1.
I took the data that started this ticket  jw03168004001_04103_00001_mirimage and processed assign_wcs using jwst version 1.16.1, 1.17.1, JWST develop Main, and JWST Develop with the branch fix using V2, V3. I did this for several programs (a science program I am involved in 4192 and a help  desk ticket that had issues with s_region (program 3224) all of them show the same thing. The s_region seemed to be fine in version 1.16.1, got messed up in 1.17.1 and then this fix using v2-v3 brings it back to version 1.16.1.  Looking at the change log I find changes to 1.17.0 (listed below that may have led to this incorrect s_region). We have a fix now for MIRI LRS fixed slit - but I do worry that one of those changes might affect other issues that we have yet to discover. 

assign_wcs (image2, spec2)


  • Use the range of points in the TabularModel to adjust the bounding_box in a

  MIRI LRS FIXEDSLIT observation.

  Ignore the bounding_box in running the inverse WCS transform when computing

  crpix. (#&#x2060;8554 <https://github.com/spacetelescope/jwst/issues/8554>_)

  • Catch NaN values in msa tables for source positions in slit and replace with

  slit center. (#&#x2060;8874 <https://github.com/spacetelescope/jwst/issues/8874>_)

  • Use pixel vertices to define s_region instead of pixel centers." (`#⁠8897

  #8897`_)

  • Fix all bounding box assignments to wcs objects so that they are correctly

  and

  specifically assigned as order F boxes. This ensures consistency with the

  assumptions made by GWCS for bounding boxes. (`#⁠8963

  #8963`_)

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Jane Morrison Huh, that's interesting.  I checkout out some of the individual SHA string commits quickly, and I think the issue must have started with #8554

The S_REGION was fine on Dec 13th, but not by Dec 19th.  Based on the logs, it looks like for years there had been a bug with NaNs in the LRS slit S_REGION, such that the pipeline aborted updating the region and kept the original (good) version from SDP:

stpipe.AssignWcsStep - INFO - There are NaNs in s_region, S_REGION not updated.

Once the NaNs in the S_REGION got fixed the region was able to be computed (incorrectly) and overwrite the correct version that was in the header.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I guess that means that the issue we've been working here must be unrelated to the original problem report, since that was filed back in September...  I'll take another look and see if there's a secondary problem.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Ok, there's two bugs.

The first bug is the issue that we've been working with assign_wcs updating the S_REGION that gets populated in the cal.fits file.  This bug is new, introduced in mid-December 2024.

There's a second bug though; S_REGION also gets updated again by the resampling code when creating the i2d.fits file.  This goes through the compute_footprint_spectral routine and seems to have some similar problems.  This bug is present in 1.16.1 (and potentially earlier), but you'll only see it if you look at the S_REGIONS in the i2d files.

We'll clearly need to fix that one too, but in general a question for Nadia Dencheva - why are we updating the S_REGION keyword so often?  If it's in RA/DEC after assign_wcs why does it need to be recreated after resampling?  (Perhaps I'm just not considering some detail about a different instrument mode).

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law Oh that makes sense - that s_region had NANs and it was not getting updated by the pipeline. We did not see the error how the S-region for MIRI LRS was determined until this was fixed. I have been going round and round on this issue because none of my MIRI LRS data had bad s_regions and that is because the s_region was likely still what set_telescope set. That was one reason I kept pointing at how assign_wcs worked because I did not see that it was a problem. Little did I know that those updates did not get applied because of Nans. Goodness that was a hidden error.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Nadia Dencheva on JIRA:

While working on the inverse transforms in #8554 I noticed that the S_REGION for MIRI LRS slit is always NaNs because the bounding box is larger than the domain of the Tabular model. The Tabular model was filling everything outside its domain with NaNs. I adjusted the bounding box to be the domain of the Tabular model. Now that the corners in V2V3 are used to compute S_REGION that change can be disregarded. 

As for why resample recomputes the S_REGION, I'm not sure. It's probably unnecessary since the WCS of the original models is not changed.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law Nadia Dencheva Melanie Clarke 
I standup today I mentioned that S_REGION is updated in resample. Melanie Clarke mentioned one reason this might be is if  the user has updated the WCS. Should we keep s_region updating for this case ? It probably does not happened much. But if we remove it - that is a sure way someone will then want do it.
Does anyone know if resampling in calspec3 would cause the S_REGION of MIRI LRS data to be larger (encompassing more than 1 slit ? 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

That's a good point; imaging data would go through tweakreg before resample for instance which certainly could update the WCS.  Dithering might be relevant here too- if multiple dithered images are combined in the resample step then the slit footprint on the sky would presumably be an amalgam of the two input slit footprints.  Whether that amalgamated footprint is meaningful or not could vary, but it should be easy enough to compute, starting from the same V2V3 corner coordinates.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

David Law Nadia Dencheva 
After resample_spec has created the drizzled data it calls the module in assign_wcs.util = update_s_region_spectral - Which is the problem for the LRS s2d s_region.

To deal with s_region in assign_wcs being incorrect I wrote a new module update_s_region_lrs - that uses the specwcs reference file.
That is not really going to work for resample. 
Could we read in the input s_regions and just use those values to find the total s_region covering the entire set of data ?
If it is just 1 file - as it would be if called from calspec2 then there is not change in s_region
If several files are used (calspec3) it would use all the s_region values to find output s_region.
Would that work ?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I think that would be consistent with what resample is doing for imaging, where it looks like the S_REGION coming out of resample for a set of dithered data is the rectangle bounding the combined S_REGION of all of the individual frames.  We'd compute the LRS slit S_REGION of the individual frames using their embedded WCS, and then figure out what new S_REGION we need to encompass all of those individual S_REGIONS.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

Nadia Dencheva For resample_spec I need to update how the final s_region is determined for MIRI LRS Fixed slit  data. Currently the same routine that caused problems in assign_wcs for MIRI LRS is being called (update_s_region_spectral).
Now resample spec seems to correctly combine the data together so all I need to do is just change how the final s_region is determined. 
Do we have a simple routine that finds the combined s_region from a set of s_regions ?
All I can find is what resample does - which is to call resample_utils.make_output_wcs
I think I could call that routine - but only use the  final s_region. It seems a bit overkill and I was wondering if you had a better idea.

 

@stscijgbot-jp
Copy link
Collaborator Author

stscijgbot-jp commented Feb 5, 2025

Comment by Jane Morrison on JIRA:

David Law Nadia Dencheva Melanie Clarke 
We have a solution for correcting the s_region when assign_wcs runs by using the V2,V3 slit corners provided by a reference file and it works well.
However correcting the s_region for the s2d files is trickier. Nadia Dencheva  suggested I try using the wcs of the output model from resample and map the bounding box x,y values to ra,dec using output wcs. That transform maps to the slit center - so we get a line. We could use the slit width (in JDOX 0.51 arc second) to fold in the width to this line. As Nadia and I were testing this new method we discoveredthat the WCS of the output of resample for calwebb_spec3 seems to be off by about a slit width. So a new problem. I have tried to summarize in two plots. In both cases I am showing the results of assign_wcs with using the V2, V3 slits corners and I plot the s_region from the cal files. In one plot I show the s_region for resample using the code on jwst/main and other plot is using the bbox of s_region, using the WCS of output and mapping to ra,dec (which CURRENTLY just maps to slit center). Look at the titles of the plots for reference.I have plotted the s_reigon for s2d from calspec2 and the final sregion from calspec3. Both show (in different ways) the sregion in calspec3 is off by about 1 slit width. For background the code on main determines the region of resample data using the min and max of ra and dec - hence the BOX for s_region.
David Law  we wanted to confirm the results that the s_region from resample calspec3 is off.

  !MIRI_LRS_slit_updates_assignwcs_no_updates_resample.png|thumbnail!
!MIRI_LRS_slit_updates.png|thumbnail!

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Jane Morrison Huh; focusing in on that second plot, it seems to be showing that the 'true' slit footprint is a good match to the slit-center line in resampled data from a single input file, but that the slit-center line when multiple images are combined together is offset from the slit-center line of the individual images.  That seems to be the case regardless of any S_REGION stuff?

Offhand I don't think that would be expected; let me try to confirm.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Ok, I think I see the issue.  Third bug?  Something is bogus with the output WCS in the s2d files, even ignoring the S_REGIONS.  I get a plot like your second one even when I just try transforming the individual XY pixels in the s2d file to RA/DEC space.  LRSwcs.png shows the slit bounding boxes for the 1st/2nd exposures (calculated in the right way from V2V3 corners and the telescope attitude in cal files), along with the RA/DEC points associated with every pixel in the individual and combined s2d files.  Clearly something is fishy- those points should match both each other and the bounding boxes.

My suspicion is that it's something to do with the rectification changing the array dimensions.

!LRSwcs.png!

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

Ok that confirms what I have found. In the PR that I have not pushed up I took the x,y corners and used the output wcs to map to ra,dec and I get the same plot you get.  am looking at the resample_spec code now. It is just so strange because the output s2d file look correct to the eye. 

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I'm poking at the resample_spec code too; so far the all_ra, all_dec and ra_center_final and dec_center_final values look ok, but something is going wrong between line 689 and 733, since the output_wcs transform is giving the wrong answers when applied to output pixel coordinates.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Maybe relevant; line 586 says

if spectral_axis == 0:  # MIRI LRS, the WCS x axis is spatial

but spectral_axis=1 when I run the code for LRS?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I don't claim to understand how the xtan and ytan projections are working, but something seems fishy with them.  If I change line 711 to

pix_to_ytan.intercept = -0.5 * (44 - 1) * pix_to_ytan.slope

(where 44 is the x size of the original array used to derive the tan plane transforms rather than the new combined array) then the across-slice position of the outputs matches almost exactly what it should.  However, the along-slice coordinates are still wrong.

@jemorrison
Copy link
Collaborator

We think alike I just did that !
It seems the wcs is set by the first data model here the entire xsize is used.
I think we are on the right track

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Christian Hayes on JIRA:

Mention of the x/ytan projections reminded me that we had an issue with this when the NIRSpec spectral resampling was updated.  I thought the methodology was copied over from MIRI, so you all may be running into the same issue we were.  We found that the resampled WCS was using the tangent plane updating either the x or y intercepts based on which one aligned better with the slit, but it needed to be updating both.  This caused an offset of the slit WCS region from where it was supposed to be

Here was the PR in case it is useful:  #8908 (and sorry if it's unrelated, but I thought I'd mention it)

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Jane Morrison on JIRA:

I have managed to get confused on axes - 0  or 1. It seems the definition is x axis = 0 and y axis =1 in code. 
At least that is how  bbox set up

 bbox  ModelBoundingBox(

    intervals={

        x0: Interval(lower=302.5, upper=346.5)

        x1: Interval(lower=6.699999999999989, upper=393.7)

    }
So in the MIRI LRS case the spectral axis is 1.

David Law not sure if you have looked in detail at xstop, x_indx, ystop y_indx 
it seems both of these are defined in terms of x_tan_array and y_tan_array which are both ~ 44 or the dimension 0 of bbox.

Does this seem correct to you ? I think so because it is spatial x and y not what x and y represent in slit.
Just wanted to double that with you

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Jane Morrison My guess is that the line that says 

if spectral_axis == 0:  # MIRI LRS, the WCS x axis is spatial

must just have an incorrect comment pasted from earlier (where it's used for if spatial_axis == 0)

I haven't looked in detail at any of the stop indices.  It looks like a number of the changes Christian Hayes mentioned for NIRSpec though affect very similar parameters, so I do wonder whether it's a similar problem.

@jemorrison jemorrison linked a pull request Feb 14, 2025 that will close this issue
10 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants