assign_wcs introduces a tilt to a line of constant wavelength

[stscijgbot-jp]

Issue JP-3491 was created on JIRA by Greg Sloan:

The current LRS distortion solutions keep a line of constant wavelength parallel to the row axis on the detector, but assign_wcs tilts the wavelengths in detector coordinates, such that the wavelength changes from one column to the next along a given row. That leads to erroneous offsets in wavelength depending on position in the cross-dispersion (i.e., spatial) direction in the spectral image. It also triggers additional bugs, including pixels with NaN values and DQ flags set to 0 instead of 1 (use instead of do_not_use), which in turn triggers other bugs. For one known dependent bug that has now been repaired, see JIRA ticket JP-3473.

[stscijgbot-jp]

Comment by David Law on JIRA:

Greg Sloan FYI relevant code starts around line 290 of https://github.com/spacetelescope/jwst/blob/master/jwst/assign_wcs/miri.py

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

Greg Sloan Please let us know the name of a suitable example dataset, for us to use when debugging.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

Any LRS slit dataset will do, as the issue is universal (same is probably true for slitless, but not confirmed). As I understand it, the existing assign_wcs takes the following steps for LRS data: 1) Fit a line to the curved trace data in the specwcs reference file (I usually refer to this file as the wavesamp file). 2) Assume a curve of constant wavelength is a line, and that it is perpendicular to the line just fitted to the curved trace data. 3) Determine the rotation necessary to align the fitted line (Step 1) to the x or y axis of the image (depending on the dispersion direction). That will rotate a line of constant wavelength to one of those axes. (I think this step is in assign_wcs, but that needs to be confirmed.)

The fundamental problem with what assign_wcs is doing is that it is ignoring most of the information in the wavesamp file. In the steps above: 1) The trace really does have curvature, and we need to preserve that and account for it. 2a) The structure of the wavesamp file also imposes the assumption that a curve of constant wavelength is a line, so that assumption in assign_wcs if fine. 2b) However, a line of constant wavelength is not necessarily perpendicular to the spectral trace at any given wavelength. In fact the current wavesamp file for the LRS was written explicitly to keep a line of constant wavelength parallel to the row axis of the detector at all wavelengths, so definitely not perpendicular to the trace, at any wavelength. 3) After the rotation, a line of constant wavelength is now tilted wrt the row axis of the detector, despite the intentional design of the row axis.

As I understand it, the rotation was implemented as a quick and necessary fix for code that did not yet exist. Unfortunately, it is not treating the data correctly. I do not have an algorithm in place that would set up an analytical transformation for the resample_spec step.

My sol'n would be numerical and would be implemented in resample_spec: First determine delta-x as a f'n of y from the wavesamp file (using x_ctr and y_ctr) and apply that shift. Then determine delta-y as a f'n of x from the wavesamp file (using the (x0, y0), (y1, y2), (etc) pairs to determine the slope of a line of constant wavelength at each x, and apply that shift.

An alternative would be to construct a grid of (x,y) pairs to define the positions of new resampled pixels on the old grid (probably almost as described above) and then use some weighted averaging scheme to determine the signal values for each new pixel (notice the waving of hands here!). That's still a numerical sol'n.

My (limited) understanding is that the assign_wcs step wants an analytical sol'n so that the software can go back and forth between detector space (x,y) and the resampled space.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

I've checked the wavelength map provided with cal.fits files (extension 4), and it shows a tilt in wavelength, so this tilt happens before resample_spec. I will post some old IDL code that we used to generate position and wavelength maps for the Spitzer/IRS wavesamp files, which have the same format as the LRS distortion files. I don't know if this code will help, as it's a workable sol'n to the problem but far from an ideal sol'n. Just thought it would be good to have it handy.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

I should offer this code with many apologies, as it's ancient and has a lot of stuff hardcoded that shouldn't be. If anyone wants to test it, I added one of the IRS wavesamp files.

[stscijgbot-jp]

Comment by David Law on JIRA:

It's way upstream of that; even miricoord gives a tilted WCS trace for the LRS as this is how the code and the corresponding assign_wcs were originally written.  This is probably 7 years ago at this point, but as I recall both miricoord and the pipeline take the same approach:

• Use the specwcs reference file to get centroid trace information
• Fit the centroid trace with a polynomial (tilted line) and use the coefficients to determine the rotation angle
• Use that rotation angle to compute the de-rotated coordinates of each pixel
• Use the derotated coordinates in spatial/spectral computations

It sounds like the main thing that needs to change is that wavelengths should be based on the actual Y detector coordinates, not the derotated Y detector coordinates.  That may be a relatively simply change.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

The spectral trace should be tilted. In fact, it should be curved, right up to when we apply resample_spec. So that's OK.

I think the thing to do is remove the rotation.

But we will need a new algorithm to be used by resample_spec.

The algorithm supplied simply produces maps of pos'n and wavelength for each pixel. It definitely doesn't mess with anything WCS.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

Nadia Dencheva  Greg Sloan I am getting up to speed on the LRS, but for other modes I have worked on we define the transformations from 'detector' to 'world'  the assign_wcs. We don't usually have an algorithm in the drizzle step or resample step that defines how to do this transformation. It seems like we need to remove the rotation in the assign_wcs but also define how to transform the detector coordinate system to the world. Nadia Dencheva is the expert on JWST transforms.  Greg could we take your algorithm which produces maps of positions and wavelengths for each pixel and convert that to a set of equations. We had to do something exactly like that for the MRS  LONG LONG ago the MRS reference file had a table of wavelengths for each pixel. We changed that to  equations. 

[stscijgbot-jp]

Comment by David Law on JIRA:

Jane Morrison I think I see how to do this; when I overhauled the LRS WCS code years back we were using a rotation to deal with the tilted trace, but it may be simpler to just use a set of shifts instead.  Experimenting with that now, crossed fingers may be able to send something to test later today.

[jemorrison]

I removed the rotation in the assign_wcs step and now I do not get the tilts in wavelength. @drlaw1558 I will see what you get after your testing.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

1) To answer Jane's question, yes, feel free to use the algorithm provided. As I recall, it has an iterative element to it, so it's not ideal.

2) To comment on David's suggestion. What he suggests is the simpler of two algorithms I was considering. Basically, shift in x in each row, then shift in y in each column. That's been my standard rectification scheme for decades. A potentially better algorithm would be to lay down some grid points for the new coordinate system on top of the old and use a sampling scheme to load data in the new grid system. Flux conservation is a concern in both cases.

(1) and (2) serve two different purposes. (1) makes maps of pos'n and wavelength for use in the original coordinate system. That would be ideal for extracting directly from the cal.fits files, and it's exactly what we would need to implement PSF-based extraction at that stage. (2) is what is needed for resample_spec.

My confusion here is why these two needs seem to be intertwined, and why we're seeing some form of resampling for the wavelengths assigned to each pixel in the cal.fits image.

[stscijgbot-jp]

Comment by David Law on JIRA:

The trick is that we needed the rotation in order to get the positions in the along-slit direction correct.  I.e., to keep the trace centroid at the same RA/DEC at all wavelengths.  It sounds like perhaps we didn't need the rotation for the wavelength axis though, if the dispersion and cross dispersion axes are not orthogonal to each other.

I've tried rewriting it to remove the rotation and instead replace it with a Y-dependent shift in the X position trace instead.  See how https://github.com/drlaw1558/jwst/tree/jp3491_lrswcs compares.

[stscijgbot-jp]

Comment by David Law on JIRA:

Just made one further update.  Apparently the pathloss step was relying on how the transforms were previous written to pull out information about the nominal reference point; I've added a dummy transform back in just so that pathloss can continue to get this information as expected.  There are definitely some differences in the s2d files with this change- not obvious to me offhand whether that's a good or a bad thing.

Just to clarify for Greg Sloan : there's no resampling of pixels at the wavelength calibration stage.  What assign_wcs needs to do is to figure out what the v2,v3,lambda is for a given detector pixel.  It does that by transforming the detector coordinates to 'slit-like' coordinates.  I.e., if you took out the spectral axis, where on the slit would you be at the reference wavelength?  Once you have that information you can use the standard imager distortion transform to compute v2/v3.

At present the mapping from actual detector pixel to 'slit' location is done by rotating the pixel coordinates by the angle of the trace, and then taking the rotated X pixel coordinate paired with the Y slit reference location.  In the modification linked above it instead does it by determining the X offset relative to the trace centroid at a given Y which will have the effect of allowing the cross-dispersion direction to be along real detector row rather than rotated rows.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

No worries, David Law, I wasn't suggesting that we were resampling the data at this stage, just the wavelength (and possibly the position) information.

This scheme of using the distortion sol'n to generate a Delta-X for each Y works for me. At some point, I will introduce a tilt to the wavelengths. Not only will it not be perpendicular to the trace, but it is likely to vary as a function of row. The sol'n would be to introduce a Delta-Y for each X and apply that shift after the Delta-X shifts. Will we be ready for that?

[stscijgbot-jp]

Comment by David Law on JIRA:

Not sure what the best means of implementing that would be, but I don't see why we couldn't adjust to it when necessary.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

David Law Greg Sloan I am still getting up to speed on this issue. As I understand it the specwcs reference file that Greg delivered in Dec included any tilt there may be by the x0...y3 regions.Combine this with the fact that we now do not treat the dispersion and along slit axis as orthogonal. David's new branch does not perform a rotation in the wavelength. But instead accounts for the fact that the spectral trace is not straight by setting up the  wavelength model using a x shift based on the specwcs xcen and ycent.

Seems reasonable. 

What I do not fully understand is how does this impact the release of the pipeline we just made. It seems the old code was either wrong or tied to a different setup of the specwcs. If it was just plain wrong we should get these change into the Pipeline ASAP. If the old specwcs was ok should we revert the specwcs until we can get these changes into pipeline ????

[stscijgbot-jp]

Comment by David Law on JIRA:

Jane Morrison My read of the situation is that what the pipeline has been doing for the past few years isn't optimal, but the latest specwcs reference file delivery didn't fundamentally change anything.  It fixed a few overall issues with the calibration data, but the x0...y3 etc information (which was in the old specwcs files too) isn't being used at all at the moment.

As I see it we've got a couple of options:

1) Remove the tilt from the wavelength solution, which is plain incorrect and has been with us for years.  This is pretty easy, and my PR above already does this.  It would set iso-lambda to along-rows, which is presumably better than the existing tilt.

2) Revise the LRS WCS code to actually use the x0...y3 information to specify the correct iso-lambda angle everywhere, possibly as some kind of 2d Tabular interpolation model.

Ideally we'd do (2).  Looking at the latest specwcs reference file though, it looks like the x0...y3 vectors currently contain no information about any iso-lambda tilts at all.  I.e., they suggest that iso-lambda is exactly along detector rows everywhere, which would functionally be identical to (1) above.  See figure attached showing a few bounding-boxes for some iso-lambda elements plotted atop the trace centroid.

I'd thus be inclined to (1) until such time as we have a reference file with non-zero tilt to actually test out.

[stscijgbot-jp]

Comment by David Law on JIRA:

(Figure of course wouldn't show tiny tilts, but numerically they're zero)

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

Yes, the tilt introduces errors in the wavelength map that I'd like to fix, but they're small compared to our overall uncertainties over most of the spectral image. They become comparable to our uncertainties at the top (i.e., at shorter wavelengths). I just checked the wavelength map for our most recent pipeline output, and here's some examples of delta lambda from one side of the spectral image to the other (i.e., col 303 to 346).

Row lambda (um) del lam 7 14.02 3.2 nm 200 10.80 4.8 nm 300 8.30 7.1 nm 350 6.49 10.1 nm 380 4.86 16.5 nm

Delta lambda between the two nods would be about 2/3 the above values, and they should average out in the Level 3 s2d.fits image.

Our uncertainty is generally +/- 20 nm at this stage, and it's probably more like 40 nm at the short-wavelength end of the spectral range.

So, yes, we need to fix this, but it's not urgent. Best to target the next feasible build.

I agree with David; option (1) is the best for the short-term, although I'm not certain what impact this will have downstream. We do need to implement option (2) long-term, because we know that we do have a small tilt that we are not yet accounting for.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

So should we merge David's PR (after a review)  into the pipeline and then try and get this new code in a release  if we need to do a patch on the current pipeline release. If we don't do a patch - we just live with it for this release ?

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

I'll defer to David here. My instinct is that this change is not important enough to justify a patch. But I'm deferring!

[stscijgbot-jp]

Comment by David Law on JIRA:

My feeling is that this doesn't necessitate a patch release; it's a small improvement to an issue that's been with us for years, so 11.0 would be fine.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

Ok sounds good. 

David do you want to make a PR using the new Branch you made. I can review it. Did you want to do more testing ? I can do any testing you think is needed

 

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

David Law I would like to move this issue forward - David do you plan to submit a PR for your changes using the above branch.

I have taken your code and I put it in a notebook and plotted what it is doing. It looks good me to. I am also using it in an LRS science program and the results are good. 

 

 

 

[stscijgbot-jp]

Comment by David Law on JIRA:

Done, see #8411

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

Fixed by https://github.com/spacetelescope/jwst/pull/8411/

[stscijgbot-jp]

Comment by David Law on JIRA:

Greg Sloan Assigning this issue to you to check out and see if the fix works as expected for you.  If you're happy enough with that already, please close the ticket.

[stscijgbot-jp]

Comment by David Law on JIRA:

Greg Sloan Bumping for you to please test if you're happy with this, ideally in the next week or so.

[stscijgbot-jp]

Comment by Ian Wong on JIRA:

David Law  I've tested the development build on several LRS slit and slitless observations, and the output isolambda lines correspond to rows. So it looks good.

[stscijgbot-jp]

Comment by Ian Wong on JIRA:

One question I had was whether or not we expect the isolambda lines for slitless mode to line up with rows in the calints images. Because no rectification is done and the spectral trace is curved, do we really expect/want the rows to have the same wavelength?

[stscijgbot-jp]

Comment by David Law on JIRA:

Thanks Ian Wong , closing.  Regarding the isolambda, my understanding from the above was that these are close-enough to lining up along rows as to not matter, even though the spectral trace is curved.  (I.e., the angle between isolambda and iso-alpha (along-slit position) is not 90 degrees and can change slightly).  I defer to the LRS team on whether or not that's accurate though.

[stscijgbot-jp]

Comment by Greg Sloan on JIRA:

The wavesamp file which produces the specwcs reference file is set to force a line of constant wavelength along the row axis of the detector, so Ian's test confirms that the assign_wcs step now preserves that geometry.

Excellent!

In answer to the questions about what really happens, yes, a line of constant wavelength is tilted slightly wrt to the row axis, and the angle of that tilt changes with row. The angle is small, and for the moment, we are fine without including it, although at some point that wil need to change. Not yet, though.