BUG: to_fits bounding_box does not work with units

[pllim]

I was looking at how to use https://gwcs.readthedocs.io/en/latest/api/gwcs.wcs.WCS.html#gwcs.wcs.WCS.to_fits . It complains that I have to give bounding box but there is no example on how to actually call that method. @mcara (and the docstring) said it is ((xmin, xmax), (ymin, ymax))

But plain floats didn't work, and passing in Quantity didn't work either.

Adapted from your example:

import gwcs
import numpy as np
from astropy import units as u
from astropy.coordinates import ICRS
from astropy.modeling import models
from gwcs import coordinate_frames as cf

shift_by_crpix = models.Shift(-(4096 - 1) * u.pix) & models.Shift(-(2048 - 1) * u.pix)
matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
                   [5.0226382102765E-06, -1.2644844123757E-05]])
rotation = models.AffineTransformation2D(matrix * u.deg, translation=[0, 0] * u.deg)
rotation.input_units_equivalencies = {"x": u.pixel_scale(1 * (u.deg / u.pix)),
                                      "y": u.pixel_scale(1 * (u.deg / u.pix))}
rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,
                                                 translation=[0, 0] * u.pix)
rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1 * (u.pix / u.deg)),
                                              "y": u.pixel_scale(1 * (u.pix / u.deg))}
tan = models.Pix2Sky_TAN()
celestial_rotation = models.RotateNative2Celestial(
    3.581704851882 * u.deg, -30.39197867265 * u.deg, 180 * u.deg)
det2sky = shift_by_crpix | rotation | tan | celestial_rotation
det2sky.name = "linear_transform"
detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix))
sky_frame = cf.CelestialFrame(reference_frame=ICRS(), name='icrs', unit=(u.deg, u.deg))
pipeline = [(detector_frame, det2sky), (sky_frame, None)]
w2 = gwcs.WCS(pipeline)  # bounding_box is None

>>> w2.to_fits(bounding_box=((0, 100), (0, 100)))
...
UnitsError: Shift: Units of input 'x', (dimensionless), could not be converted to required input units of pix (unknown)

>>> w2.to_fits(bounding_box=([0, 100] * u.pix, [0, 100] * u.pix))
...
UnitConversionError: Can only apply 'less' function to dimensionless quantities when other argument
is not a quantity (unless the latter is all zero/infinity/nan)

[pllim]

@mcara said the reported errors is a bug, so I updated the issue title.

[nden]

Does it work if you use a bounding_box define with units (in px)?

[pllim]

As stated above, I tried w2.to_fits(bounding_box=([0, 100] * u.pix, [0, 100] * u.pix)) and it still errored out. Did I do it wrong?

[nden]

bounding_box=([[0, 100], [0, 100]] * u.pix) may be better but I may be misremembering. @WilliamJamieson can confirm if this is a bug.

[WilliamJamieson]

bounding_box=([[0, 100], [0, 100]] * u.pix) may be better but I may be misremembering. @WilliamJamieson can confirm if this is a bug.

This doesn't work either. It produces this error:

TypeError: type Quantity doesn't define __round__ method

Traceback:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 w2.to_fits(bounding_box=bounding_box)

File ~/projects/gwcs/gwcs/gwcs/wcs.py:2415, in WCS.to_fits(self, bounding_box, max_pix_error, degree, max_inv_pix_error, inv_degree, npoints, crpix, projection, bin_ext_name, coord_col_name, sampling, verbose)
   2412         degree = 1
   2413         max_inv_pix_error = None
-> 2415     hdr = self._to_fits_sip(
   2416         celestial_group=celestial_group,
   2417         keep_axis_position=True,
   2418         bounding_box=bounding_box,
   2419         max_pix_error=max_pix_error,
   2420         degree=degree,
   2421         max_inv_pix_error=max_inv_pix_error,
   2422         inv_degree=inv_degree,
   2423         npoints=npoints,
   2424         crpix=crpix,
   2425         projection=projection,
   2426         matrix_type='PC-CDELT1',
   2427         verbose=verbose
   2428     )
   2429     use_cd = 'A_ORDER' in hdr
   2431 else:

File ~/projects/gwcs/gwcs/gwcs/wcs.py:1787, in WCS._to_fits_sip(self, celestial_group, keep_axis_position, bounding_box, max_pix_error, degree, max_inv_pix_error, inv_degree, npoints, crpix, projection, matrix_type, verbose)
   1785 # 0-based crpix:
   1786 if crpix is None:
-> 1787     crpix1 = round(bb_center[input_axes[0]], 1)
   1788     crpix2 = round(bb_center[input_axes[1]], 1)
   1789 else:

[pllim]

Still a problem. Very frustrating.

Maybe you should remove units from your Getting Started example, since a lot of stuff (distortion, this) won't work with the units.

[mairanteodoro]

@pllim What worked for me was to use dimensionless units for the shift and affine transformations:

import gwcs
import numpy as np
from astropy import units as u
from astropy.coordinates import ICRS
from astropy.modeling import models
from gwcs import coordinate_frames as cf

shift_by_crpix = models.Shift(-(4096 - 1)) & models.Shift(-(2048 - 1))  # <-- removed units from here
matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
                   [5.0226382102765E-06, -1.2644844123757E-05]])
rotation = models.AffineTransformation2D(matrix, translation=[0, 0])  # <-- removed units from here
rotation.input_units_equivalencies = {"x": u.pixel_scale(1 * (u.deg / u.pix)),
                                      "y": u.pixel_scale(1 * (u.deg / u.pix))}
rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,
                                                 translation=[0, 0] * u.pix)
rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1 * (u.pix / u.deg)),
                                              "y": u.pixel_scale(1 * (u.pix / u.deg))}
tan = models.Pix2Sky_TAN()
celestial_rotation = models.RotateNative2Celestial(
    3.581704851882 * u.deg, -30.39197867265 * u.deg, 180 * u.deg)
det2sky = shift_by_crpix | rotation | tan | celestial_rotation
det2sky.name = "linear_transform"
detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix))
sky_frame = cf.CelestialFrame(reference_frame=ICRS(), name='icrs', unit=(u.deg, u.deg))
pipeline = [(detector_frame, det2sky), (sky_frame, None)]
w2 = gwcs.WCS(pipeline)  # bounding_box is None

After those changes, running w2.to_fits(bounding_box=((0, 100), (0, 100))) will result in

(NAXIS   =                    2 / number of array dimensions                     
 NAXIS1  =                  101                                                  
 NAXIS2  =                  101                                                  
 WCSAXES =                    2 / Number of coordinate axes                      
 WCSNAME = 'icrs    '                                                            
 CRPIX1  =                 51.0 / Pixel coordinate of reference point            
 CRPIX2  =                 51.0 / Pixel coordinate of reference point            
 PC1_1   = 1.29087916909098E-05 / Coordinate transformation matrix element       
 PC1_2   = 5.94419672962468E-06 / Coordinate transformation matrix element       
 PC2_1   = 5.01416888582657E-06 / Coordinate transformation matrix element       
 PC2_2   = -1.2648737798687E-05 / Coordinate transformation matrix element       
 CDELT1  =                  1.0 / Coordinate increment at reference point  
 CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  
 CUNIT1  = 'deg'                / Units of coordinate increment and value        
 CUNIT2  = 'deg'                / Units of coordinate increment and value        
 CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
 CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
 CRVAL1  =        3.50740873226 / [deg] Coordinate value at reference point      
 CRVAL2  =     -30.387022471349 / [deg] Coordinate value at reference point      
 LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
 LATPOLE =     -30.387022471349 / [deg] Native latitude of celestial pole        
 MJDREF  =                  0.0 / [d] MJD of fiducial time                       
 RADESYS = 'ICRS'               / Equatorial coordinate system                   
 SIPMXERR= 1.38956422009025E-06 / Max diff from GWCS (equiv pix).                
 COMMENT FITS WCS created by approximating a gWCS                                ,
 [])