Add common resample code to stcal. #320

mcara · 2024-11-27T16:54:03Z

Closes #279

This PR add the common resample code used by both JWST and Roman pipelines to stcal. Also, for the first time, this PR adopts the new drizzle API from spacetelescope/drizzle#134 for the resample code used in the pipelines. This PR is a second attempt to #279 and therefore merging this PR should close #279.

This work is related to https://jira.stsci.edu/browse/AL-835

Tasks

update or add relevant tests
update relevant docstrings and / or docs/ page
Does this PR change any API used downstream? (if not, label with no-changelog-entry-needed)
- write news fragment(s) in changes/: echo "changed something" > changes/<PR#>.<changetype>.rst (see below for change types)
- run regression tests with this branch installed ("git+https://github.com/<fork>/stcal@<branch>")
  - jwst regression test
  - romancal regression test

news fragment change types...

changes/<PR#>.apichange.rst: change to public API
changes/<PR#>.bugfix.rst: fixes an issue
changes/<PR#>.general.rst: infrastructure or miscellaneous change

codecov · 2024-11-27T16:57:04Z

Codecov Report

Attention: Patch coverage is 54.79323% with 481 lines in your changes missing coverage. Please review.

Project coverage is 83.59%. Comparing base (35726e9) to head (c0b9971).
Report is 5 commits behind head on main.

Files with missing lines	Patch %	Lines
src/stcal/resample/resample.py	60.90%	355 Missing ⚠️
src/stcal/resample/utils.py	16.55%	126 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main     #320       +/-   ##
===========================================
+ Coverage   29.57%   83.59%   +54.01%     
===========================================
  Files          36       52       +16     
  Lines        7949    10076     +2127     
===========================================
+ Hits         2351     8423     +6072     
+ Misses       5598     1653     -3945

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

braingram

I left a few comments to try and address some of the CI failures.

There are also check-type failures.

I also don't see any tests. Are there tests that can be moved here from jwst/romancal?

src/stcal/resample/resample.py

kmacdonald-stsci

Overall it looks fine, but it's being reported that a lot of the code isn't tested by any test. Maybe add tests to get more code coverage.

src/stcal/resample/resample.py

kmacdonald-stsci · 2024-11-27T23:51:49Z

src/stcal/resample/resample.py

+        attributes = {
+            "data",
+            "wcs",
+            "wht",


Is this short for "weight"? When I see "wht", I think "what". I recommend simply using "weight".

These have direct analogue in the JWST DataModel. That is why I preserved the name.

Here are a couple of examples:

https://github.com/spacetelescope/romancal/blob/e250097fb7113be808836a796dc73ba7f687c7b8/romancal/resample/gwcs_drizzle.py#L78

https://github.com/spacetelescope/jwst/blob/40e43c8dbe5c3cec269c4d014dfbf53ddb23a3bb/jwst/resample/resample.py#L413

https://github.com/spacetelescope/jwst/blob/363a9ad90eb66111bcecd65edab36b2a9e631feb/jwst/resample/resample.py#L499

It looks like this is being changed in jwst, so can it also be changed here?

src/stcal/resample/resample.py

emolter

I agree with Brett, the major high-level thing I see is that this needs unit tests

src/stcal/resample/__init__.py

emolter · 2024-11-29T14:53:02Z

src/stcal/resample/resample.py

+                    int(output_wcs.bounding_box[0][1] + halfpix),
+                )
+            else:
+                # TODO: In principle, we could compute footprints of all


Can this TODO (or any other ones in this code) be removed and if desired, turned into a GitHub issue or JIRA ticket so we can track work still needed?

It could, but does it make sense to open a ticket to modify code that does not yet exist in the repository?

It's not necessarily my decision but I would say yes, because this is definitely going to get added at some point and it's easier to track issues than code comments

src/stcal/resample/utils.py

emolter

looks good overall, most of my comments are surface-level

src/stcal/resample/resample.py

emolter · 2024-12-04T23:20:57Z

src/stcal/resample/resample.py

+            'pixmap': pixmap,
+            'scale': iscale,
+            'weight_map': weight,
+            'wht_scale': 1.0,  # hard-coded for JWST count-rate data


Is this okay for Romancal? If not, perhaps now is the time to allow it as a variable

seems to be OK: https://github.com/spacetelescope/romancal/blob/e250097fb7113be808836a796dc73ba7f687c7b8/romancal/resample/gwcs_drizzle.py#L206

great, thanks for checking. maybe update the comment to say that it was also tested for Roman, otherwise it looks like it may have been accidentally ported over from jwst without additional checking

maybe update the comment to say that it was also tested for Roman, otherwise it looks like it may have been accidentally ported over from jwst without additional checking

I did not test it and I am not aware of this being tested but I assume the resample code was tested as a whole by Roman. I actually do not see under what circumstances wht_scale would need to be != 1. I think the comment should be removed as I think units of the data are irrelevant here.

src/stcal/resample/resample.py

perrygreenfield

Most of the comments regard utils.py

src/stcal/resample/utils.py

perrygreenfield · 2024-12-10T23:28:57Z

src/stcal/resample/utils.py

+
+
+def compute_wcs_pixel_area(wcs, shape=None):
+    """ Computes pixel area in steradians.


I'm not happy with this and the following routine. It is I think overly complex and doesn't even explain the methodology of how the pixel area is computed. I gather that it is an average area of pixels over the whole valid image area. But why not simply compute the area of one pixel in the center? If that isn't useful, then it needs to explain why it must be the average. So at the very least this is very poorly explained, in my opinion. Perhaps it can be handled in a subsequent PR if this is urgent. At the very least the docstring should be more descriptive for this PR.

The purpose of this PR is to move code that may be shared by multiple pipelines (JWST, Roman, ...) to stcal. This PR attempts to preserve how current pipelines work and I do not think this should be modified here.

I'm not happy with this and the following routine. It is I think overly complex and doesn't even explain the methodology of how the pixel area is computed.

I will add a little more documentation.

But why not simply compute the area of one pixel in the center?

We would not need to compute anything if we knew where/how PIXAR_* were computed and them maybe we could used them directly without computing anything.

I gather that it is an average area of pixels over the whole valid image area. But why not simply compute the area of one pixel in the center?

It is too complicated to describe it here but this came from the ambiguity of how PIXAR_SR was defined for JWST images. I attempted to provide a longer explanation here: spacetelescope/jwst#7894 (comment) and maybe this: spacetelescope/jwst#7894 (comment)

Further reading: spacetelescope/jwst#7668 (comment) and other comments in that issue.

Perhaps it can be handled in a subsequent PR if this is urgent.

Again, this is code already used in JWST. Please open an issue if you feel this is incorrect.

src/stcal/resample/utils.py

perrygreenfield

LGTM

schlafly

I left a few minor comments inline, but this looks good, thanks Mihai.

schlafly · 2024-12-16T17:09:26Z

src/stcal/resample/resample.py

+            # TODO: if we want to support adding more data to
+            # existing output models, we need to also store weights
+            # for variance arrays:
+            # var_rnoise_weight
+            # var_flat_weight
+            # var_poisson_weight


Do you just need one of these (i.e., for squared kernel terms), or do you need all of them?

schlafly · 2024-12-16T17:20:55Z

src/stcal/resample/utils.py

+        xmax = min(data_shape[1] - 1, int(x2 + 0.5))
+        ymax = min(data_shape[0] - 1, int(y2 + 0.5))
+
+    return xmin, xmax, ymin, ymax


Feel free to ignore, but I'll note that I'm not immediately sure what to make of the usage of whole pixels vs. half pixels. I think I interpret that the bbox has the corners of the pixels, with values like (-0.5, N - 0.5). In that case the mins at int(x + 0.5) make sense, but the maxes at int(x + 0.5) look a bit weird. I might have guessed floor(x + 0.5), ceil(x - 0.5). Of course, there's not a lot of space between int(x + 0.5) and ceil(x - 0.5); maybe just a comment about what the intent is re half / whole pixels.

Somehow, github won't allow me to reply separately to each comment.

Do you just need one of these (i.e., for squared kernel terms), or do you need all of them?

I think so. I do not think it needs to be done now, especially that I've heard that Roman is no longer interested in this functionality, but this is something to keep in mind for later only if we will keep current method of estimating variances. If we propagate errors like it was proposed in spacetelescope/jwst#8997 and spacetelescope/drizzle#163, then these arrays are not needed.

the maxes at int(x + 0.5) look a bit weird. I might have guessed floor(x + 0.5), ceil(x - 0.5).

int(x + 0.5) is the same as int(floor(x + 0.5)). xmin, xmax, ymin, ymax must all be integer numbers. If we make a rounding error here, it preferably should be to expand the box, not to shrink it: minor efficiency loss is probably preferable to missing data.

src/stcal/resample/utils.py

pyproject.toml

Co-authored-by: Ned Molter <[email protected]>

emolter

Mostly grammar, but a few other things I'm wondering about too.

Also, we should figure out why the oldestdeps test is running into astropy version conflicts and failing to build.

emolter · 2025-01-31T16:07:13Z

docs/conf.py

+nitpicky = False
+
+nitpick_ignore = [
+    ('py:obj', 'type'),
+]


What required this change?

emolter · 2025-01-31T16:09:14Z

docs/stcal/resample/description.rst

+sum.  The process is:
+
+#. For each input model, take the square root of each of the read noise variance
+   array to make an error image.


Suggested change

array to make an error image.

arrays to make an error image.

emolter · 2025-01-31T16:09:48Z

docs/stcal/resample/description.rst

+
+#. Square the resampled read noise to make a variance array.
+
+   a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight


Suggested change

a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight

#. If the resampling `weight_type` is an inverse variance map (`ivm`), weight

emolter · 2025-01-31T16:11:38Z

docs/stcal/resample/description.rst

+in the ``con`` attribute in the output data model . Each pixel in the context
+image is a bit field that encodes
+information about which input image has contributed to the corresponding
+pixel in the resampled data array. Context image uses 32 bit integers to encode


Suggested change

pixel in the resampled data array. Context image uses 32 bit integers to encode

pixel in the resampled data array. The context image uses 32 bit integers to encode

emolter · 2025-01-31T16:11:47Z

docs/stcal/resample/description.rst

+image is a bit field that encodes
+information about which input image has contributed to the corresponding
+pixel in the resampled data array. Context image uses 32 bit integers to encode
+this information and hence it can keep track of only 32 input images.


Suggested change

this information and hence it can keep track of only 32 input images.

this information, and hence it can keep track of only 32 input images.

emolter · 2025-01-31T16:18:30Z

src/stcal/resample/resample.py

+       scale factors needed to obtain correctly convert resampled counts to
+       fluxes.
+    3. For each input image computes coordinate transformations (``pixmap``)
+       from coordinate system of the input image to the coordinate system of


Suggested change

from coordinate system of the input image to the coordinate system of

from the coordinate system of the input image to the coordinate system of

emolter · 2025-01-31T16:18:54Z

src/stcal/resample/resample.py

+    3. For each input image computes coordinate transformations (``pixmap``)
+       from coordinate system of the input image to the coordinate system of
+       the output image.
+    4. For each input image computes weight image.


Suggested change

4. For each input image computes weight image.

4. Computes the weight image for each input image.

emolter · 2025-01-31T16:32:25Z

src/stcal/resample/resample.py

+    def get_input_model_pixel_area(self, model):
+        """
+        Computes or retrieves pixel area of an input model. Currently,
+        this is the average pixel area of input model's pixels within either


Suggested change

this is the average pixel area of input model's pixels within either

this is the average pixel area of the input model's pixels within either

emolter · 2025-01-31T16:33:34Z

src/stcal/resample/resample.py

+        difference in the definition of the pixel area reported in model's
+        ``meta`` and the pixel area at the location used to construct


Suggested change

difference in the definition of the pixel area reported in model's

``meta`` and the pixel area at the location used to construct

difference in the definition of the pixel area reported in

``model.meta`` and the pixel area at the location used to construct

...although, now that this is written in such a way that everything is passed in and out with dictionaries, is meta still relevant?

emolter · 2025-01-31T16:33:47Z

src/stcal/resample/resample.py

+        ``meta`` and the pixel area at the location used to construct
+        output WCS from the WCS of input models using ``pixel_scale_ratio``.
+
+        Intensity scale factor is computed elsewhere as the ratio of the value


Suggested change

Intensity scale factor is computed elsewhere as the ratio of the value

The intensity scale factor is computed elsewhere as the ratio of the value

mcara requested a review from a team as a code owner November 27, 2024 16:54

mcara self-assigned this Nov 27, 2024

mcara requested review from braingram, perrygreenfield, nden and emolter November 27, 2024 16:54

mcara added the resample label Nov 27, 2024

mcara mentioned this pull request Nov 27, 2024

Move common resample code to stcal spacetelescope/jwst#8986

Open

10 tasks

mcara force-pushed the resample-common-code2 branch from 50484f2 to ebbc003 Compare November 27, 2024 17:22

braingram reviewed Nov 27, 2024

View reviewed changes

src/stcal/resample/resample.py Outdated Show resolved Hide resolved

src/stcal/resample/resample.py Outdated Show resolved Hide resolved

src/stcal/resample/resample.py Outdated Show resolved Hide resolved

kmacdonald-stsci approved these changes Nov 28, 2024

View reviewed changes

emolter reviewed Nov 29, 2024

View reviewed changes

mcara force-pushed the resample-common-code2 branch from 119a118 to d20777f Compare December 3, 2024 18:31

github-actions bot added the installation label Dec 3, 2024

emolter reviewed Dec 4, 2024

View reviewed changes

src/stcal/resample/utils.py Show resolved Hide resolved

emolter reviewed Dec 4, 2024

View reviewed changes

perrygreenfield reviewed Dec 13, 2024

View reviewed changes

github-actions bot added the documentation Improvements or additions to documentation label Dec 15, 2024

mcara force-pushed the resample-common-code2 branch 2 times, most recently from 1c4b44c to 6743256 Compare December 16, 2024 07:07

perrygreenfield approved these changes Dec 16, 2024

View reviewed changes

schlafly approved these changes Dec 16, 2024

View reviewed changes

braingram mentioned this pull request Jan 27, 2025

Hackday mteodoro sky variance readnoise spacetelescope/romancal#1556

Open

10 tasks

braingram reviewed Jan 27, 2025

View reviewed changes

pyproject.toml Outdated Show resolved Hide resolved

mcara added 5 commits January 28, 2025 01:17

Add common resample code to stcal.

7c6fb53

remove unused imports

7cd2ad8

remove available memory check

2822226

Changes to accomodate outlier detection and spec in jwst

4659a8d

make code more robust

08d5f1e

mcara and others added 8 commits January 28, 2025 01:17

fix type

8406f03

fix imports

eef3952

address reviewers' comments - part 1

d557826

Update src/stcal/resample/resample.py

0f180ee

Co-authored-by: Ned Molter <[email protected]>

Update src/stcal/resample/resample.py

a6c836b

Co-authored-by: Ned Molter <[email protected]>

add more docs

4e33c2e

Fix type warning

5849a65

Use updated drizzle to deal with gwcs bbox

969ab74

mcara force-pushed the resample-common-code2 branch from 6743256 to 969ab74 Compare January 28, 2025 07:37

extra checks in resample_range

a05ebbd

emolter reviewed Jan 31, 2025

View reviewed changes

braingram mentioned this pull request Jan 31, 2025

Use astropy pixel_to_pixel in calc_pixmap spacetelescope/drizzle#161

Open

fix a bug in resample_range

c0b9971



		def compute_wcs_pixel_area(wcs, shape=None):
		""" Computes pixel area in steradians.


		#. Square the resampled read noise to make a variance array.

		a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight

	a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight
	#. If the resampling `weight_type` is an inverse variance map (`ivm`), weight

	pixel in the resampled data array. Context image uses 32 bit integers to encode
	pixel in the resampled data array. The context image uses 32 bit integers to encode

	this information and hence it can keep track of only 32 input images.
	this information, and hence it can keep track of only 32 input images.

	from coordinate system of the input image to the coordinate system of
	from the coordinate system of the input image to the coordinate system of

	4. For each input image computes weight image.
	4. Computes the weight image for each input image.

	this is the average pixel area of input model's pixels within either
	this is the average pixel area of the input model's pixels within either

		difference in the definition of the pixel area reported in model's
		``meta`` and the pixel area at the location used to construct

	Intensity scale factor is computed elsewhere as the ratio of the value
	The intensity scale factor is computed elsewhere as the ratio of the value

Add common resample code to stcal. #320

Are you sure you want to change the base?

Add common resample code to stcal. #320

Conversation

mcara commented Nov 27, 2024

Tasks

codecov bot commented Nov 27, 2024 • edited Loading

Codecov Report

braingram left a comment

Choose a reason for hiding this comment

kmacdonald-stsci left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emolter left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mcara Dec 13, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emolter left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

perrygreenfield left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mcara Dec 15, 2024 • edited Loading

Choose a reason for hiding this comment

perrygreenfield left a comment

Choose a reason for hiding this comment

schlafly left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

emolter left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Nov 27, 2024 •

edited

Loading

mcara Dec 13, 2024 •

edited

Loading

mcara Dec 15, 2024 •

edited

Loading

emolter left a comment •

edited

Loading