astropy/astropy — pull request #1960 by mdboom (2014-01-13 18:30:44+00:00)
Title: 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.