segmentation fault for astropy.wcs.Wcsprm.sub()

Hi, 

I am facing a problem that I can not solve with astropy (v3), so here it is. I am reading a large fits file produced by a radio astronomical software (AIPS++) with a NAXIS4 image, of which two dimensions are collapsed.

astropy.io.fits read the image and header properly, using:
> data, hdr =  astropy.io.fits.getdata(filename, 0, header=True)

, I then create a WCS object with
> wcs = astropy.wcs.WCS(hdr)

The following is the result of wcs.to_header() : 

```
WCSAXES =                    4 / Number of coordinate axes                      
CRPIX1  =                 8193 / Pixel coordinate of reference point            
CRPIX2  =                 8193 / Pixel coordinate of reference point            
CRPIX3  =                    1 / Pixel coordinate of reference point            
CRPIX4  =                    1 / Pixel coordinate of reference point            
CDELT1  =   -0.000555555555556 / [deg] Coordinate increment at reference point  
CDELT2  =    0.000555555555556 / [deg] Coordinate increment at reference point  
CDELT3  =                    1 / Coordinate increment at reference point        
CDELT4  =             33333332 / [Hz] Coordinate increment at reference point   
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CUNIT4  = 'Hz'                 / Units of coordinate increment and value        
CTYPE1  = 'RA---SIN'           / Right ascension, orthographic/synthesis project
CTYPE2  = 'DEC--SIN'           / Declination, orthographic/synthesis projection 
CTYPE3  = 'STOKES'             / Coordinate type code                           
CTYPE4  = 'FREQ'               / Frequency (linear)                             
CRVAL1  =        35.3324166667 / [deg] Coordinate value at reference point      
CRVAL2  =       -4.51685555556 / [deg] Coordinate value at reference point      
CRVAL3  =                    1 / Coordinate value at reference point            
CRVAL4  =        322601561.836 / [Hz] Coordinate value at reference point       
LONPOLE =                  180 / [deg] Native longitude of celestial pole       
LATPOLE =       -4.51685555556 / [deg] Native latitude of celestial pole        
RESTFRQ =            306000000 / [Hz] Line rest frequency                       
RESTWAV =                    0 / [Hz] Line rest wavelength                      
EQUINOX =                 2000 / [yr] Equinox of equatorial coordinates         
SPECSYS = 'TOPOCENT'           / Reference frame of spectral coordinates        
MJD-OBS =                55794 / [d] MJD of observation matching DATE-OBS       
DATE-OBS= '2011-08-21T00:00:00.000000' / ISO-8601 observation date matching MJD-
```

I wanted to "cleverly" remove the extra empty dimensions (3 & 4) by doing : 
> goodAxis = [axis['coordinate_type'] == 'celestial' for axis in wcs.get_axis_types()]
> subData   = np.squeeze(data, axis=np.where(np.array(goodAxis) == True))
and updating the header with 
> subWCS   = wcs.sub([astropy.wcs.WCSSUB_CELESTIAL])

however if I force a garbage collection ( by gc.collect() after the wcs.sub() ) I have a nice segmentation fault (python 2.7, ipython v0.13.2)

Any idea ?

a.