lyalpha / spalipy Goto Github PK

This is for older versions of astropy. For astropy<=3, we have the code and error:

>>> from spalipy import spalipy
>>> s= spalipy.Spalipy("r0204859_UT4-median.cat", "r0204805_UT4-median.cat", "r0204859_UT4-median.fits", output_filename='r0204859_UT4-median-aligned.fits')
>>> s.main()
Matched 8429 detections within 5 pixels with affine transformation
Affine alignment pixel residuals [median (stddev)]: x = 0.047 (0.402), y = -0.035 (0.461)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.5/dist-packages/spalipy/spalipy.py", line 180, in main
    self.align()
  File "/usr/local/lib/python3.5/dist-packages/spalipy/spalipy.py", line 394, in align
    self.source_fits.writeto(output_filename,
AttributeError: 'list' object has no attribute 'writeto'

where r0204859_UT4-median.fits should be read in as an HDUList but is read in instead as a "regular" list.

I think a simple fix would be to add the lines

if type(self.source_fits) is list:
            self.source_fits = fits.HDUList(self.source_fits)

at 393 prior to

self.source_fits.writeto(output_filename,
                                     overwrite=overwrite)

In some cases the psf is so poor that SExtractor splits a star into multiple objects (e.g. if it's a donut or signicantly extended/curved). This causes the splines to sometimes go crazy. Although not tested, this is probably going to occur in extremely crowded fields, or if there are extreme seeing ratios too.

[Ideally SExtractor should just be run with correct parameters to stop this happening, but corner cases occur.]

One solution is to take only the nearest matched detection between the source (after affine transformation) and template detections only. But this would result in a biased transform.
This could be alleviated by maybe choosing the brightest detection within the match distance, but again this may change per detection match.

The detections with multiple matches within the matching distance could be thrown out (a la alipy) but this is at the cost of tie points for the fits.

Need to look into how the matching from alipy is actually done. Can source detections 'share' a match in the template or must they be unique? Could it be possible to combine the multiple matches to a single, flux weighted match based on the X_IMAGE, Y_IMAGE, FLUX_BEST of the matches? Likely a speed hit if this route pursued.

Trying to align r0034872_UT2-median.fits to r0033604_UT2-median.fits is a good example of this problem.

As a related aside - the splines shifts should be clipped to the min/max of the coordinate residuals since there shouldn't be any correction larger than this amount.

There may be nan or inf entries lingering in the detection table(s) after sep.extract() calls. I think this is probably the cause of some warnings like:

RuntimeWarning: invalid value encountered in greater_equal result = getattr(super(), op)(other)

that have been reported to arise.

Simply pruning any rows where any(np.is_finite(row)) would be enough. This is a bit annoying to implement in astropy.Table (cf. pandas.DataFrame.drop_na() for example) as one might have to loop over rows or columns.
(Note we want to remove the rows, not just replace the non-finite values)

Instead of requiring SExtractor catalogues to be passed, the code should fall-back to a simple source detection routine (e.g. photutils) and do its own extraction.

Splines are prone to going wild in the absence of constraints and/or conflicted cross matches (potentially caused by 6f705cd removing all catalogue detections in a badly distorted corner of the image, for example).

An additional argument function should be passable to Spalipy that takes spline or poly as a value and implements the correct function.

• spline_order -> order and can be used by both
• references to spline_transform etc. change to warp_transform

Astropy 2DPolynomial model seem most appropriate and can be easily fitted. The fitting can be a least squares in x and y, mimicking the spline procedure.

Use the remaining coordinate residuals (after affine+spline transformation) to determine some form of 'goodness' statistic for the transformation.

Add a results method that provides feedback on:

• shift, scale, rotation
• min/median/max spline shifts in x and y
• rms of the transformed coordinate residual offsets
  With option to plot:
• Images of the x and y spline surfaces
• An image of the residual offsets to see the spatial dependence.

Since scipy.interpolation.affine_transform fails when any nan or inf are present, options could be added to handle them.

A simple replacement_val argument can be used to replace them:

source_data[~np.isfinite] = replacement_val

(An error message should be returned if replacement_val is not set and np.sum(np.isnan(source_data)) > 0.)

More complexly replacement_val could be something like "median", "bilinear", "bicubic", which would denote an interpolation method to replace them.

Alternatively have a nan_mask option, nan_mask == ~np.isfinite(source_data), which has the transformations applied to it and returned to or written for the user. Any pixel >0 indicates a pixel affected by a nan (or inf).

Add logger (or at least print statements) to provide feedback on how the run is going.

Combine the affine transformation with the spline shifts in order to only have one stage of interpolation.