[AstroPy] Rebin FITS images and preserve or recalculate WCS?

Tue Aug 4 10:04:12 EDT 2020

Hi,

There is a rebin method in MPDAF, with an option for the margins, and 
using the same reshape trick as Gary: Image.rebin, which calls the 
generic ND method on the parent class, and WCS.rebin (class based on 
Astropy one):
https://github.com/musevlt/mpdaf/blob/master/lib/mpdaf/obj/image.py#L2380
https://github.com/musevlt/mpdaf/blob/master/lib/mpdaf/obj/data.py#L1377
https://github.com/musevlt/mpdaf/blob/master/lib/mpdaf/obj/coords.py#L1579

Simon


On 8/4/20 9:37 AM, Bernstein, Gary M wrote:
> Hi Eric - I actually just faced the rebinning issue with a summer student and had produced this function, which is very general about the dimensions of the input array.  It also handles the not-an-even-multiple case by trimming the input image to be a multiple of the binning factor.  The overflow issue is a good point; if your inputs are ushort then you might need to convert the type, which could be done at the “copy()” line below.
> 
> And I see that the SIP polynomials are applied to the pixel coordinates rather than the intermediate world coordinates, which I didn’t know about, so indeed that complicates the WCS a bit.  Looking at https://irsa.ipac.caltech.edu/data/SPITZER/docs/files/spitzer/shupeADASS.pdf
> it seems like the polynomial arguments are still relative to CRPIXn, so the CRPIX transformation in the earlier message would still work and the CDp_q scaling still needed.  But you’d have to alter all the [AB]_p_q keywords by a factor s**(p+q-1) and the inverse [AB]P_p_q by s**(1-p-q).  So it’s still not too messy, but you do have to go hunting through the header for things now, & I understand your desired to find a previously vetted solution!
> 
> Gary
> 
> 
> def blockavg(array,factor):
>      '''Average the array in blocks of `factor` pixels on each axis.
>      `factor` can be an integer, in which case it is applied to
>      all axes, or it can be an array, in which case each axis
>      is reduced by corresponding element of factor until either
>      axes or factor array are exhausted.
>      
>      If array dimension is not divisible by factor, it will
>      be clipped to make it so.
>      
>      Returns the reduced array.
>      '''
>      
>      if type(factor)==int:
>          # If factor is an int, apply to all dims
>          ff = [factor]*array.ndim
>      else:
>          # Just copy each element of an iterable
>          ff = [i for i in factor]
>          
>      naxes = min(array.ndim,len(ff))  # Number of axes to avg
> 
>      # First I'll truncate all axes to be divisible by factors:
>      s = [ slice((array.shape[i]//ff[i])*ff[i]) for i in range(naxes)]
>      tmp = array[tuple(s)].copy()
> 
>      # Now reduce each axis in turn
>      for i in range(naxes):
>          s = tmp.shape[:i] + (tmp.shape[i]//ff[i],ff[i]) + tmp.shape[i+1:]
>          tmp = np.sum( tmp.reshape(s),axis=i+1)
>      
>      # Divide to turn into average - delete this line if you want to sum...
>      tmp = tmp / np.prod(ff)
>      return tmp
> 
> 
>> On Aug 4, 2020, at 9:16 AM, Eric Jensen <ejensen1 at swarthmore.edu> wrote:
>>
>> Thanks, Gary and Petr!  I appreciate the answers so far.  A couple of additional notes:
>>
>> I should have been clear that I realize that the basics of the binning part are hard with numpy operations, though I appreciate the clear code from Gary.   There are possibly edge-of-image issues, but those are probably best handled by just restricting the binning factor to numbers that divide evenly into the array dimensions.
>>
>> However, there are potential bit-depth issues - if I add together several pixels with 50,000 counts (previously stored as 16-bit unsigned integers) I’ll need to change the bit depth.  Maybe astropy’s FITS-handling routines make that automatic when writing an array to a new FITS file.
>>
>> That said, it’s mostly the subtleties of the WCS part that I’d like to handle with vetted tools if possible.
>>
>> For example, if I look at a FITS header from an image solved by astrometry.net, it has not just CRVAL1,2 + CRPIX1,2 + a CD matrix, but it also has SIP distortion coefficients.  It’s possible that those are negligible for our images (I haven’t yet checked them to work out the details of the standard to figure out how much difference they make) and they could just be ignored.  (And our images straight off the telescope, using MaximDL, have some kind of proprietary/undocumented distortion correction scheme that I’ll just need to ignore anyway, I think, so maybe the distortion issue isn’t important.  I’ll have to test that.)
>>
>> But it highlights the general issue that WCS has lots of permutations, and it seems better to handle it (if possible) with tested tools such as already exist in astropy.  The question is whether there’s a useful way to deploy the astropy (or other library) WCS functionality here, or whether it is indeed best just to do it manually with the linear coefficients following the solution outlined by Gary and Petr already.
>>
>> One specific question:
>>
>>> On Aug 4, 2020, at 6:33 AM, Petr Kubánek <petr at kubanek.net> wrote:
>>>
>>> And if you happen to have WCS in CD_ matrix, divide it by 2 (scalar, eg. all members).
>>
>> These would need to be *multiplied* by the binning factor, yes?
>>
>> If there is no general solution already implemented and tested, I certainly can dig into it more myself - I think the basic outline (minus distortion) is already given, modulo subtleties I haven’t though of.  But in general I’ve encouraged my students to use tested libraries rather than re-inventing the wheel when possible, so I’m trying to take my own advice!  :-)
>>
>> Thanks for any additional thoughts,
>>
>> Eric
>>
>>
>>> On Aug 4, 2020, at 6:33 AM, Petr Kubánek <petr at kubanek.net> wrote:
>>>
>>> And if you happen to have WCS in CD_ matrix, divide it by 2 (scalar, eg. all members).
>>>
>>> Petr
>>>
>>>> 4. 8. 2020 v 11:44 dop., Bernstein, Gary M <garyb at PHYSICS.UPENN.EDU>:
>>>>
>>>> Hi Eric -
>>>> I don’t know whether this exists but in this case it might be faster to write than to do the Google search.  The block summing can be done in numpy by reshaping the array.  In the case where the original image shape is divisible by the shrinkage factor s, you’d do
>>>>
>>>> m,n = img.shape
>>>> img = np.sum( img.reshape(m//s,s,n), axis=1)  # Contract along y
>>>> Img = np.sum( img.reshape(m//s,n//s,s),axis=2). # Contract along x
>>>>
>>>> And to fix the WCS, you need to divide CRPIXn by s, and multiply CDELTn by s.  There is some slight complication to the first operation because FITS assumes 1-indexed pixels, and the integer marks the center of a pixel, so its really CRPIXn -> (CRPIXn—0.5)/s + 0.5
>>>>
>>>> I think this should work universally because the linear transform is always the first step of a WCS mapping.
>>>>
>>>> Cheers,
>>>> Gary
>>>>
>>>>> On Aug 3, 2020, at 11:01 PM, Eric Jensen <ejensen1 at swarthmore.edu> wrote:
>>>>>
>>>>> Hi all,
>>>>>
>>>>> We are looking at purchasing a new CMOS camera that has 4-micron pixels, which would significantly undersample our PSF. (The camera is otherwise excellent for our needs, e.g. excellent QE.)  So I’m looking at whether there would be a straightforward way in scripted post-processing to reduce the image file sizes while not losing information.  The existing camera driver doesn’t support more than 2x2 binning at this time.
>>>>>
>>>>> Is there available Python code that can take a FITS image with a valid WCS in the header and rebin it, and output a FITS image that still has a valid WCS for the rebinned image?  The rebinning should also preserve flux, though I think that’s easier than the WCS part.  Nothing obvious turns up in a little searching, but I could easily have missed it.
>>>>>
>>>>> It seems that Montage might be able to do this (mShrink) but if possible I’d prefer a pure-Python solution, as it would be easier to implement under Windows on our observatory control computer.  (A fallback would be to install a Linux distro on top of the Windows Subsystem for Linux, but if it’s possible to do it purely in Python that would be a lot simpler.)
>>>>>
>>>>> Thanks for any thoughts,
>>>>>
>>>>> Eric
>>>>>
>>>>> Eric Jensen
>>>>> Professor of Astronomy
>>>>> Dept. of Physics and Astronomy
>>>>> Swarthmore College
>>>>> _______________________________________________
>>>>> AstroPy mailing list
>>>>> AstroPy at python.org
>>>>> https://mail.python.org/mailman/listinfo/astropy
>>>>
>>>> _______________________________________________
>>>> AstroPy mailing list
>>>> AstroPy at python.org
>>>> https://mail.python.org/mailman/listinfo/astropy
>>>
>>> _______________________________________________
>>> AstroPy mailing list
>>> AstroPy at python.org
>>> https://mail.python.org/mailman/listinfo/astropy
>>
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at python.org
>> https://mail.python.org/mailman/listinfo/astropy
> 
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy
>