Currently Warp will only takes an input image and produces an output pixel grid. However, astronomical data are often accompanied with their standard deviation image (containing the error in each pixel that has been kept/propagated from the very start of the reduction).

Within a pipeline, this "error" image also needs to be warped with the input (and then stacked) so we can propagate the errors of each observation to the final stacked image.

In the current implementation (as of version 0.19), a user will have to manually create a variance image from their error image (because errors get added after being taken to the square root of 2), then run Warp on it (similar to how they warp their input).

However, the WCS calculation of the pixel edges will be the same between the input and error images, so will the polygon clipping and fractional area measurement. So it is a waste of time, RAM and CPU to run Warp twice in such scenarios.

It is more efficient for Warp to also have an option called '--instd' and '--stdhdu' (similar to MakeCatalog) to let the user provide an STD image. When the user provides this image (it is optional), Warp will read that image, take it to the power of 2 (to become variance), and then do the warping also on the variance image in the exact same way it works on the input's pixel values. It will finally take the square root of the final variance in each output pixel and return the error (standard deviation) in each output pixel as a separate HDU in the output file.

Of course, the power-of-2 and the square-root operations can take extra CPU, so like MakeCatalog, we can have a '--variance' option that would lead it to assume that the error image is already the variance, and it would also just output the variance of the output (not standard deviation).

This should be pretty easy to implement, and we can test it with the 'SKY_STD' extension of NoiseChisel's outputs on the image that are to be warped (like the example in the book).