Error in ResolutionRenderer for observations with different WCS

[eramey16]

Hello,
I am having an issue making a scarlet frame where the input images have different WCS. I've been working with Charlotte Ward on this (using the time-domain branch), but so far we haven't been able to solve it, so I figured I'd make sure it was documented as an issue. One potential cause we identified was that the WCS of my images seems to be in a different format than in your test cases (using CD rather than PC values), although just converting the WCS to use PC parameters doesn't solve the problem. I've put together a script that demonstrates the error outside of our pipeline that I can send by email, but I'll include a summary below.

Here is how I am making clips of our larger image files:

# Open image file
    with fits.open(img_file, memmap=True) as file:
        img_data = file[0].data
        img_hdr = file[0].header
        img_wcs = WCS(file[0].header)
        x, y = img_wcs.world_to_pixel(sc) # sc is the coordinates of the transient
    # Make a clip of the image
    med = np.median(img_data)
    clip = Cutout2D(img_data, sc, size=101, wcs=img_wcs, mode='partial', fill_value=med)

I do the same to the weight files, except the fill value is zero for edge cases. The PSFs are loaded from PSFex files (all are different (square) shapes, but even if I truncate them to the same shape I get the same behavior). Then, later, I load the data into the scarlet observations like this:

# Create the scarlet1 observation object and add it to our list of observations
    wcs = clip.wcs
    obs = scarlet.Observation(data,
        wcs=wcs,
        psf=psf,
        channels=channel,
        weights=weight)
    obssingle.append(obs)

scarlet.display.show_observation works fine, and all clips seem to be aligned to within a pixel. Finally, here is the code that produces the error:

model_psf_s = scarlet.GaussianPSF(sigma=0.9)
model_frame_s = scarlet.Frame.from_observations([obssingle[0], obssingle[1]], 
                                                coverage='intersection',
                                                model_psf=model_psf_s)
print(model_frame_s)

And the error itself:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[14], line 2
      1 model_psf_s = scarlet.GaussianPSF(sigma=0.9)
----> 2 model_frame_s = scarlet.Frame.from_observations([obssingle[0], obssingle[1]], 
      3                                                 coverage='intersection',
      4                                                 model_psf=model_psf_s)
      5 print(model_frame_s)

File /usr/local/lib/python3.10/dist-packages/scarlet-1.0.1+g84193b5-py3.10-linux-x86_64.egg/scarlet/frame.py:291, in Frame.from_observations(observations, model_psf, model_wcs, obs_id, coverage)
    289 # Match observations to this frame
    290 for obs in observations: 
--> 291     obs.match(model_frame)
    293 return model_frame

File /usr/local/lib/python3.10/dist-packages/scarlet-1.0.1+g84193b5-py3.10-linux-x86_64.egg/scarlet/observation.py:109, in Observation.match(self, model_frame, renderer)
    105                 self.renderer = ConvolutionRenderer(
    106                     self, model_frame, convolution_type="fft"
    107                 )
    108             else:
--> 109                 self.renderer = ResolutionRenderer(self, model_frame)
    110 else:
    111     assert isinstance(renderer, Renderer)

File /usr/local/lib/python3.10/dist-packages/scarlet-1.0.1+g84193b5-py3.10-linux-x86_64.egg/scarlet/renderer.py:354, in ResolutionRenderer.__init__(self, data_frame, model_frame, padding)
    351     axes = (1, 2)
    353 # Computes the resampling/convolution matrix
--> 354 resconv_op = self.sinc_shift(self.diff_kernel, self.shifts, axes)
    355 self._resconv_op = np.array(resconv_op, dtype=model_frame.dtype) * self.h ** 2
    357 if self.isrot:

File /usr/local/lib/python3.10/dist-packages/scarlet-1.0.1+g84193b5-py3.10-linux-x86_64.egg/scarlet/renderer.py:444, in ResolutionRenderer.sinc_shift(self, imgs, shifts, axes)
    442     shifter = np.array(interpolation.mk_shifter(self._fft_shape, real=True))
    443 else:
--> 444     shifter = np.array(interpolation.mk_shifter(self._fft_shape))
    445 # Shift
    446 if 0 in axes:
    447     # Fourier shift

File /usr/local/lib/python3.10/dist-packages/autograd/numpy/numpy_wrapper.py:58, in array(A, *args, **kwargs)
     56 t = builtins.type(A)
     57 if t in (list, tuple):
---> 58     return array_from_args(args, kwargs, *map(array, A))
     59 else:
     60     return _array_from_scalar_or_array(args, kwargs, A)

File /usr/local/lib/python3.10/dist-packages/autograd/tracer.py:48, in primitive.<locals>.f_wrapped(*args, **kwargs)
     46     return new_box(ans, trace, node)
     47 else:
---> 48     return f_raw(*args, **kwargs)

File /usr/local/lib/python3.10/dist-packages/autograd/numpy/numpy_wrapper.py:77, in array_from_args(array_args, array_kwargs, *args)
     75 @primitive
     76 def array_from_args(array_args, array_kwargs, *args):
---> 77     return _np.array(args, *array_args, **array_kwargs)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.

From digging into this a bit, it seems that the build_diffkernel method in frame.py returns a non-square array, although I can't say why this might be happening, as all the clips and PSFs are square. I also noticed that the get_angles function in interpolation.py returns one element of the list as its own array, but I can't find evidence that this is causing the main issue.

Thanks for your help!
-- Emily