Reprojection performance

[wtbarnes]

Currently, I'm using reproject_interp to do the reprojection form the spectral cube to the overlappogram. This is exceedingly slow and even seems to run out of memory when run on my laptop. Presumably this is due to the sheer number of arguments passed to map_coordinates which seems to be the bottleneck for doing this type of computation. There are a few possible strategies for speeding this up

• embarrasing parallelism--reproject each wavelength individually and then reassemble
• use dask-image to parallelize map_coordinates call
• use cupy to offload map_coordinates computation to the GPU

[wtbarnes]

It seems that the performance bottleneck is actually in the pixel-to-pixel transformation and not in the map_coordinates step, at least in the cases I've looked at thus far.

My solution for now is to do each wavelength slice separately and to attempt to do this in parallel. This makes each pixel-to-pixel transformation far less onerous and means we can pass it off to something like Dask to make this calculation more scalable.

There is currently a draft of this in the sandbox.ipynb folder. I would like to add this is as a separate block in the reproject_to_overlappogram function in the overlappy repo so that applying this operation in parallel is just a matter of setting a flag (and instantiating a Dask client).

[wtbarnes]

Some of the work being done on reproject could also be helpful here: https://github.com/astropy/reproject/pull/376

[wtbarnes]

And here: https://github.com/astropy/reproject/pull/374

[astrofrog]

Happy to discuss this use case more as I'd be interested in making sure it works fast! Just to check, are you reprojecting a 3D spectral cube just in the spatial dimension or are you doing a full 3D reprojection changing both celestial and spectral WCS?

[wtbarnes]

Great! In general, the use case here is the latter: a full 3D reprojection in both the spatial and wavelength dimensions.

[astrofrog]

But just to check, are the celestial and spectral WCSes decoupled in your case?

[wtbarnes]

One is decoupled and one is not. The use case here is reprojecting from a space-space-wavelength cube to a WCS where the wavelength dimension is coupled to one (or sometimes both) of the spatial dimensions.

In general, the PCij for the decoupled spectral cube is,

| cos(alpha) -sin(alpha) 0 |
| sin(alpha) cos(alpha)  0 |
|     0           0      1 |

while the coupled "overlappogram" WCS is

| cos(alpha) -sin(alpha) -b*cos(alpha-gamma) |
| sin(alpha) cos(alpha)  -b*sin(alpha-gamma)|
|     0           0          1 |

where alpha is the roll angle, gamma is the angle of the dispersive grating, and b is the spectral order.

For a simple case of alpha=gamma=0, b=1, this just reduces to

| 1 0 -1 |
| 0 1  0 |
| 0 0  1 |

a skew between the first spatial dimension.

I'm on mobile right now, but tomorrow I can provide full example headers for both WCSs.