Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

#102 Add multi-image FITS I/O functions #144

Merged
merged 11 commits into from

4 participants

@rmjarvis
Owner

Four new functions and some mild clean up of the previous FITS I/O functions:

  • These read/write multi-extension fits files with each image in a separate hdu:
readMulti(fits)
writeMulti(image_list, fits, add_wcs=True, clobber=True)
  • These read/write fits data cubes. i.e. a 3-d array all in one hdu:
readCube(fits)
writeCube(image_list, fits, add_wcs=True, clobber=True):
  • Removed the confusing read as a static method of Image and ImageView classes. This was confusing because it doesn't necessarily return the type specified. Rather the data type is determined from the fits file being read, not the class being used to call the function. Everyone who used the function used the galsim.fits.read syntax, so it was pretty painless to remove.
  • Added the CD matrix to the WCS being written to the fits headers. Considering that this was extolled as an advantage to keeping the pixel size in the Image classes, it seemed like we should actually implement it.
  • Added unit tests for the write function as well as all the new ones. (Only read had been unit tested, comparing the read values to a reference image.)
@rmjarvis
Owner

Not sure why it putting the read and write functions on the same line. Obviously those are two different functions in each case.

@barnabytprowe

Hi Mike - am I right in interpreting the changes here that the external FITS files used for testing have been removed, and instead the FITS images are internally generated?

If so, I have concerns that this damages the generality of the unit tests. Those FITS images were made externally (IDL) and this fact alone was enough to catch an important issue (the little-endian / big-endian stuff, resulting in the case catching you see in fits.py). There may be more we haven't come across yet - I don't know, but I'm not prepared to bet against it.

So, if anything, I'd prefer to see more independently-generated FITS test images in the repo for these tests, from alternative sources/systems, rather than fewer. They are tiny files.

@rmjarvis
Owner

That's fine. The problem with the old ones was they were square. So that misses any errors in the nx vs ny, which are easy to make. I had one in my readCube function, for example. So if you want to make fits files from IDL that match these, I'm fine with adding that unit test back in.

@rmjarvis
Owner

I guess I should mention that I don't have IDL nor even know how to use it. Otherwise I'd do that myself.

@rmandelb
Owner

I agree with Barney's concern about unit tests. But otherwise, it looks good to me.

@barnabytprowe

OK Mike, fair enough, I'll make fits files you need and push them.

@barnabytprowe

Hi all, I've put the new IDL-generated FITS files into the repo.

I don't want to put the IDL code used to make them into the repo, not least as it has other dependencies (IDL Astrolib). So, just in case I want to do this again, I'm going to take the liberty of repeating it here for posterity and in all it's clunky IDL glory...

PRO make_test_ims

; Single images
test_array = [[11, 21, 31, 41, 51, 61, 71], $
              [12, 22, 32, 42, 52, 62, 72], $
              [13, 23, 33, 43, 53, 63, 73], $
              [14, 24, 34, 44, 54, 64, 74], $ 
              [15, 25, 35, 45, 55, 65, 75]]
writefits, 'testS.fits', test_array
writefits, 'testI.fits', long(test_array)
writefits, 'testF.fits', float(test_array)
writefits, 'testD.fits', double(test_array)

; Then do cubes with NIMAGES = 12 and multi-extension FITS images with
; each extension
; having ext_no added to it
nimages = 12
test_cube = intarr([size(test_array, /DIM), 12])
for k=0, nimages-1 do begin

   test_cube[*, *, k] = test_array + k

endfor

; First write these cubes out
writefits, 'test_cubeS.fits', test_cube
writefits, 'test_cubeI.fits', long(test_cube)
writefits, 'test_cubeF.fits', float(test_cube)
writefits, 'test_cubeD.fits', double(test_cube)

; Multi-ext: start with 16 bit ints
filename = 'test_multiS.fits'
spawn, "rm "+filename
mkhdr, header, test_array, /EXTEND
writefits, filename, test_array, header
for k=1, nimages-1 do begin

   writefits, filename, test_cube[*, *, k], /APPEND

endfor
; then proceed to other types
filename = 'test_multiI.fits'
spawn, "rm "+filename
mkhdr, header, long(test_array), /EXTEND
writefits, filename, long(test_array), header
for k=1, nimages-1 do begin

   writefits, filename, long(test_cube[*, *, k]), /APPEND

endfor
;
filename = 'test_multiF.fits'
spawn, "rm "+filename
mkhdr, header, float(test_array), /EXTEND
writefits, filename, float(test_array), header
for k=1, nimages-1 do begin

   writefits, filename, float(test_cube[*, *, k]), /APPEND

endfor
;
filename = 'test_multiD.fits'
spawn, "rm "+filename
mkhdr, header, double(test_array), /EXTEND
writefits, filename, double(test_array), header
for k=1, nimages-1 do begin

   writefits, filename, double(test_cube[*, *, k]), /APPEND

endfor
END

Right, I'm very happy to merge this. There was one thing I noted: that readMulti crashes when reading in multiple extension fits where the primary HDU is empty (i.e. if all the data has gone into the extensions). I'm not sure if we ever expect data like that, but I do know that I inadvertently made some for the tests and it took a little bit of hamfisted IDL to work.

If you guys think this is worth it, I'll open an issue to think about this and maybe address it.

@rmjarvis
Owner

I think I would put this in the repo if you think it will be useful to make similar things down the line. It doesn't have to be runnable from galsim prerequisites. I would even call it a .txt file to make it clear this isn't meant to be runnable code. Just a guide for running something like it in the future.

@rmjarvis
Owner

As for the empty primary HDU thing, I don't mind adding that as a feature. We can check the first one, and if empty then skip it. (Is it testable as if not hdu in python? Or is it more complicated, like there is an hdu but the data array comes back as zero size?

@rmandelb
Owner

A few comments:

1) everything compiles, tests pass, etc. on my machine
2) I'm fine with merging
3) I agree with Mike: we could even make a directory for these extra utilities... e.g., I have an IDL script from this morning that I used to make the postage stamps and catalogs that it would be good to have here on github. I'm sure that when we start doing validation / comparison with imengine and the GREAT08 and GREAT10 codes, we will have scripts and things that have other dependencies that we might want to save as well.

@rmjarvis
Owner

I just looked at the diff for your commit. I don't like removing the tests of the write functions. These need to be unit tested as well. I thought you were just going to add unit tests with the IDL files as well. The write functions should write to a different name and check that reading the image back in works as expected.

@barnabytprowe
@barnabytprowe

(But it won't be for ~1.5 hours as I really do have to go into the office today, and the commute takes a little time.) Once that's done, are you guys happy to merge? I'll let you know.

@rmjarvis
Owner

No rush. The important bits are already merged into #103 where I needed it. And I don't think anyone else is waiting on this.

@reikonakajima
Collaborator

Everything complies and tests successfully on my Mac OS X 10.7.4 and Ubuntu 10.04.4.

@barnabytprowe

Right, changes pushed and tests run for me, now with Mike's old as well as mine based on the FITS in the repo. I added the IDL script to devutils/external/, which a brief comment about it's use as well as a ref. to this Issue.

Merge at will, as far as I'm concerned!

@rmjarvis
Owner

OK. Looks good to me now. Tests passed on Mac and linux for me too.

@rmjarvis rmjarvis merged commit 0759776 into master
@rmandelb rmandelb deleted the #102 branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
76 devutils/external/make_test_ims.pro
@@ -0,0 +1,76 @@
+PRO make_test_ims
+
+; IDL script used to generate external test images for GalSim.
+; Used by tests/test_Image.py, see Issue #144.
+;
+
+; Single images
+test_array = [[11, 21, 31, 41, 51, 61, 71], $
+ [12, 22, 32, 42, 52, 62, 72], $
+ [13, 23, 33, 43, 53, 63, 73], $
+ [14, 24, 34, 44, 54, 64, 74], $
+ [15, 25, 35, 45, 55, 65, 75]]
+writefits, 'testS.fits', test_array
+writefits, 'testI.fits', long(test_array)
+writefits, 'testF.fits', float(test_array)
+writefits, 'testD.fits', double(test_array)
+
+; Then do cubes with NIMAGES = 12 and multi-extension FITS images with
+; each extension
+; having ext_no added to it
+nimages = 12
+test_cube = intarr([size(test_array, /DIM), 12])
+for k=0, nimages-1 do begin
+
+ test_cube[*, *, k] = test_array + k
+
+endfor
+
+; First write these cubes out
+writefits, 'test_cubeS.fits', test_cube
+writefits, 'test_cubeI.fits', long(test_cube)
+writefits, 'test_cubeF.fits', float(test_cube)
+writefits, 'test_cubeD.fits', double(test_cube)
+
+; Multi-ext: start with 16 bit ints
+filename = 'test_multiS.fits'
+spawn, "rm "+filename
+mkhdr, header, test_array, /EXTEND
+writefits, filename, test_array, header
+for k=1, nimages-1 do begin
+
+ writefits, filename, test_cube[*, *, k], /APPEND
+
+endfor
+; then proceed to other types
+filename = 'test_multiI.fits'
+spawn, "rm "+filename
+mkhdr, header, long(test_array), /EXTEND
+writefits, filename, long(test_array), header
+for k=1, nimages-1 do begin
+
+ writefits, filename, long(test_cube[*, *, k]), /APPEND
+
+endfor
+;
+filename = 'test_multiF.fits'
+spawn, "rm "+filename
+mkhdr, header, float(test_array), /EXTEND
+writefits, filename, float(test_array), header
+for k=1, nimages-1 do begin
+
+ writefits, filename, float(test_cube[*, *, k]), /APPEND
+
+endfor
+;
+filename = 'test_multiD.fits'
+spawn, "rm "+filename
+mkhdr, header, double(test_array), /EXTEND
+writefits, filename, double(test_array), header
+for k=1, nimages-1 do begin
+
+ writefits, filename, double(test_cube[*, *, k]), /APPEND
+
+endfor
+END
+
View
233 galsim/fits.py
@@ -21,9 +21,9 @@ def write(image, fits, add_wcs=True, clobber=True):
If add_wcs evaluates to True, a 'LINEAR' WCS will be added using the Image's
bounding box. This is not necessary to ensure an Image can be round-tripped
through FITS, as the bounding box (and scale) are always saved in custom header
- keys. If add_true is a string, this will be used as the WCS name.
+ keys. If add_wcs is a string, this will be used as the WCS name.
- This function called be called directly as "galsim.fits.write(image, ...)",
+ This function can be called directly as "galsim.fits.write(image, ...)",
with the image as the first argument, or as an image method: "image.write(...)".
Setting clobber=True when 'fits' is a string will silently overwrite existing files.
@@ -34,11 +34,13 @@ def write(image, fits, add_wcs=True, clobber=True):
hdus = fits
else:
hdus = pyfits.HDUList()
+
if len(hdus) == 0:
hdu = pyfits.PrimaryHDU(image.array)
else:
hdu = pyfits.ImageHDU(image.array)
hdus.append(hdu)
+
hdu.header.update("GS_SCALE", image.scale, "GalSim Image scale")
hdu.header.update("GS_XMIN", image.xMin, "GalSim Image minimum X coordinate")
hdu.header.update("GS_YMIN", image.xMin, "GalSim Image minimum Y coordinate")
@@ -56,12 +58,140 @@ def write(image, fits, add_wcs=True, clobber=True):
"coordinate system value at reference pixel")
hdu.header.update("CRPIX1" + wcsname, 1, "coordinate system reference pixel")
hdu.header.update("CRPIX2" + wcsname, 1, "coordinate system reference pixel")
+ hdu.header.update("CD1_1" + wcsname, image.scale, "CD1_1 = pixel_scale")
+ hdu.header.update("CD2_2" + wcsname, image.scale, "CD2_2 = pixel_scale")
+ hdu.header.update("CD1_2" + wcsname, 0, "CD1_2 = 0")
+ hdu.header.update("CD2_1" + wcsname, 0, "CD2_1 = 0")
if isinstance(fits, basestring):
if clobber and os.path.isfile(fits):
os.remove(fits)
hdus.writeto(fits)
+def writeMulti(image_list, fits, add_wcs=True, clobber=True):
+ """
+ Write the image to a multi-extension FITS file:
+ - If 'fits' is a pyfits.HDUList, the images will be appended as new HDUs.
+ The user is responsible for calling fits.writeto(...) afterwards.
+ - If 'fits' is a string, it will be interpreted as a filename for a new
+ multi-extension FITS file.
+
+ If add_wcs evaluates to True, a 'LINEAR' WCS will be added using each Image's
+ bounding box. This is not necessary to ensure an Image can be round-tripped
+ through FITS, as the bounding box (and scale) are always saved in custom header
+ keys. If add_wcs is a string, this will be used as the WCS name.
+
+ Setting clobber=True when 'fits' is a string will silently overwrite existing files.
+ """
+ import pyfits # put this at function scope to keep pyfits optional
+
+ if isinstance(fits, pyfits.HDUList):
+ hdus = fits
+ else:
+ hdus = pyfits.HDUList()
+
+ for image in image_list:
+ write(image, hdus, add_wcs=add_wcs, clobber=clobber)
+
+ if isinstance(fits, basestring):
+ if clobber and os.path.isfile(fits):
+ os.remove(fits)
+ hdus.writeto(fits)
+
+
+def writeCube(image_list, fits, add_wcs=True, clobber=True):
+ """
+ Write the image to a FITS file as a data cube:
+ - If 'fits' is a pyfits.HDUList, the cube will be appended as new HDUs.
+ The user is responsible for calling fits.writeto(...) afterwards.
+ - If 'fits' is a string, it will be interpreted as a filename for a new
+ FITS file.
+ Normally 'image_list' is a python array of Image's (or ImageViews).
+ Each is required to have the same size (nx,ny). If not a ValueError is raised.
+
+ The image_list can also be either an array of numpy arrays or a 3d numpy array,
+ in which case this is written to the fits file directly. In the former case,
+ not explicit check is made that the numpy arrays are all the same shape, but
+ a numpy exception will be raised which we let pass upstream unmolested.
+
+ If add_wcs evaluates to True, a 'LINEAR' WCS will be added using the first Image's
+ bounding box. No check is made to confirm that all images have the
+ same origin and pixel scale. If add_wcs is a string, this will be used as the WCS name.
+
+ Setting clobber=True when 'fits' is a string will silently overwrite existing files.
+ """
+ import numpy
+ import pyfits # put this at function scope to keep pyfits optional
+
+ if isinstance(fits, pyfits.HDUList):
+ hdus = fits
+ else:
+ hdus = pyfits.HDUList()
+
+ try:
+ cube = numpy.asarray(image_list)
+ nimages = cube.shape[0]
+ nx = cube.shape[1]
+ ny = cube.shape[2]
+ # Use default values for xmin, ymin, scale
+ scale = 1
+ xMin = 1
+ yMin = 1
+ except:
+ nimages = len(image_list)
+ if (nimages == 0):
+ raise IndexError("In writeCube: image_list has no images")
+ im = image_list[0]
+ nx = im.xMax - im.xMin + 1
+ ny = im.yMax - im.yMin + 1
+ scale = im.scale
+ xMin = im.xMin
+ yMin = im.yMin
+ array_shape = (nimages, ny, nx)
+ cube = numpy.array([[[]]])
+ cube.resize(array_shape)
+ for k in range(nimages):
+ im = image_list[k]
+ nx_k = im.xMax-im.xMin+1
+ ny_k = im.yMax-im.yMin+1
+ if nx_k != nx or ny_k != ny:
+ raise IndexError("In writeCube: image %d has the wrong shape"%k +
+ "Shape is (%d,%d). Should be (%d,%d)"%(nx_k,ny_k,nx,ny))
+ cube[k,:,:] = image_list[k].array
+
+ if len(hdus) == 0:
+ hdu = pyfits.PrimaryHDU(cube)
+ else:
+ hdu = pyfits.ImageHDU(cube)
+ hdus.append(hdu)
+
+ hdu.header.update("GS_SCALE", scale, "GalSim Image scale")
+ hdu.header.update("GS_XMIN", xMin, "GalSim Image minimum X coordinate")
+ hdu.header.update("GS_YMIN", xMin, "GalSim Image minimum Y coordinate")
+
+ if add_wcs:
+ if isinstance(add_wcs, basestring):
+ wcsname = add_wcs
+ else:
+ wcsname = ""
+ hdu.header.update("CTYPE1" + wcsname, "LINEAR", "name of the coordinate axis")
+ hdu.header.update("CTYPE2" + wcsname, "LINEAR", "name of the coordinate axis")
+ hdu.header.update("CRVAL1" + wcsname, xMin,
+ "coordinate system value at reference pixel")
+ hdu.header.update("CRVAL2" + wcsname, yMin,
+ "coordinate system value at reference pixel")
+ hdu.header.update("CRPIX1" + wcsname, 1, "coordinate system reference pixel")
+ hdu.header.update("CRPIX2" + wcsname, 1, "coordinate system reference pixel")
+ hdu.header.update("CD1_1" + wcsname, scale, "CD1_1 = pixel_scale")
+ hdu.header.update("CD2_2" + wcsname, scale, "CD2_2 = pixel_scale")
+ hdu.header.update("CD1_2" + wcsname, 0, "CD1_2 = 0")
+ hdu.header.update("CD2_1" + wcsname, 0, "CD2_1 = 0")
+
+ if isinstance(fits, basestring):
+ if clobber and os.path.isfile(fits):
+ os.remove(fits)
+ pyfits.writeto(fits,cube)
+
def read(fits):
"""
@@ -78,18 +208,8 @@ def read(fits):
Not all FITS pixel types are supported (only those with C++ Image template
instantiations are: short, int, float, and double).
- This function can be called directly as "galsim.fits.read(...)", or as a static
- method of an image class: "ImageD.read(...)". Note, however, that in the
- latter case the image type returned is determined by the type of the FITS file,
- not the image class (in other words, "ImageD.read(...)" might return an ImageF).
+ This function is called as "im = galsim.fits.read(...)"
"""
- # MJ: I find this last syntax: ImageD.read(...) a bit confusing, since as you
- # point out the return value isn't necessarily the class you call it from.
- # Also, the return value is now an ImageView, not an Image, but that should
- # be transparent to the user.
- # So I'd recommend removing the ImageD.read(...) syntax and just having
- # the galsim.fits.read(...) syntax.
-
import pyfits # put this at function scope to keep pyfits optional
if isinstance(fits, basestring):
@@ -120,14 +240,95 @@ def read(fits):
image.scale = scale
return image
-# inject read/write as methods of Image classes
+def readMulti(fits):
+ """
+ Construct an array of ImageViews from a FITS data cube.
+ - If 'fits' is a pyfits.HDUList, it will read images from these
+ - If 'fits' is a string, it will be interpreted as a filename to open and read
+
+ If the FITS header has GS_* keywords, these will be used to initialize the
+ bounding box and scale. If not, the bounding box will have (xMin,yMin) at
+ (1,1) and the scale will be set to 1.0.
+
+ Not all FITS pixel types are supported (only those with C++ Image template
+ instantiations are: short, int, float, and double).
+ """
+
+ import pyfits # put this at function scope to keep pyfits optional
+
+ if isinstance(fits, basestring):
+ hdu_list = pyfits.open(fits)
+ elif isinstance(fits, pyfits.HDUList):
+ hdu_list = fits
+ else:
+ raise TypeError("In readMulti, fits is not a string or HDUList")
+
+ image_list = []
+ for hdu in hdu_list:
+ image_list.append(read(hdu))
+ return image_list
+
+def readCube(fits):
+ """
+ Construct an array of ImageViews from a FITS data cube.
+ - If 'fits' is a pyfits.HDUList, the Primary HDU will be used.
+ - If 'fits' is a pyfits.PrimaryHDU or pyfits.ImageHDU, that HDU will be used.
+ - If 'fits' is a string, it will be interpreted as a filename to open;
+ the Primary HDU of that file will be used.
+
+ If the FITS header has GS_* keywords, these will be used to initialize the
+ bounding boxes and scales. If not, the bounding boxes will have (xMin,yMin) at
+ (1,1) and the scale will be set to 1.0.
+
+ Not all FITS pixel types are supported (only those with C++ Image template
+ instantiations are: short, int, float, and double).
+
+ This function is called as "image_list = galsim.fits.readCube(...)"
+ """
+ import pyfits # put this at function scope to keep pyfits optional
+
+ if isinstance(fits, basestring):
+ fits = pyfits.open(fits)
+ if isinstance(fits, pyfits.HDUList):
+ fits = fits[0]
+
+ xMin = fits.header.get("GS_XMIN", 1)
+ yMin = fits.header.get("GS_YMIN", 1)
+ scale = fits.header.get("GS_SCALE", 1.0)
+ pixel = fits.data.dtype.type
+
+ try:
+ Class = _galsim.ImageView[pixel]
+ except KeyError:
+ raise TypeError("No C++ Image template instantiation for pixel type %s" % pixel)
+ # Check through byteorder possibilities, compare to native (used for numpy and our default) and
+ # swap if necessary so that C++ gets the correct view.
+ if fits.data.dtype.byteorder == '!':
+ if native_byteorder == '>':
+ pass
+ else:
+ fits.data.byteswap(True)
+ elif fits.data.dtype.byteorder in (native_byteorder, '=', '@'):
+ pass
+ else:
+ fits.data.byteswap(True) # Note inplace is just an arg, not a kwarg, inplace=True throws
+ # a TypeError exception in EPD Python 2.7.2
+
+ nimages = fits.data.shape[0]
+ image_list = []
+ for k in range(nimages):
+ image = Class(array=fits.data[k,:,:], xMin=xMin, yMin=yMin)
+ image.scale = scale
+ image_list.append(image)
+ return image_list
+
+
+# inject write as methods of Image classes
for Class in _galsim.Image.itervalues():
Class.write = write
- Class.read = staticmethod(read)
for Class in _galsim.ImageView.itervalues():
Class.write = write
- Class.read = staticmethod(read)
for Class in _galsim.ConstImageView.itervalues():
Class.write = write
View
1  tests/Image_comparison_images/.gitignore
@@ -0,0 +1 @@
+*_internal.fits
View
BIN  tests/Image_comparison_images/testD.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/testF.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/testI.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/testS.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_cubeD.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_cubeF.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_cubeI.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_cubeS.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_multiD.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_multiF.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_multiI.fits
Binary file not shown
View
BIN  tests/Image_comparison_images/test_multiS.fits
Binary file not shown
View
206 tests/test_Image.py
@@ -11,11 +11,11 @@
Each image is 5x7 pixels^2 and if each pixel is labelled (x, y) then each pixel value is 10*x + y.
The array thus has values:
-15 25 35 45 55
-14 24 34 44 54
-13 23 33 43 53 ^
-12 22 32 42 52 |
-11 21 31 41 51 y
+15 25 35 45 55 65 75
+14 24 34 44 54 64 74
+13 23 33 43 53 63 73 ^
+12 22 32 42 52 62 72 |
+11 21 31 41 51 61 71 y
x ->
@@ -23,6 +23,11 @@
checks, etc.
Images are in S, I, F & D flavours.
+
+There are also four FITS cubes, and four FITS multi-extension files for testing. Each is 12
+images deep, with the first image being the reference above and each subsequent being the same
+incremented by one.
+
"""
try:
@@ -47,6 +52,8 @@
[14, 24, 34, 44, 54, 64, 74],
[15, 25, 35, 45, 55, 65, 75] ]).astype(np.int16)
+nimages = 12 # Depth of FITS datacubes and multi-extension FITS files
+
datadir = os.path.join(".", "Image_comparison_images")
def test_Image_basic():
@@ -131,32 +138,175 @@ def test_Image_FITS_IO():
"""Test that all four FITS reference images are correctly read in by both PyFITS and our Image
wrappers.
"""
- # The test fits file was made with the following array. I didn't remake the file
- # when I changed the ref_array to have ncol != nrow.
- fits_ref_array = np.array([
- [00, 10, 20, 30],
- [01, 11, 21, 31],
- [02, 12, 22, 32],
- [03, 13, 23, 33]]).astype(np.int16)
for i in xrange(ntypes):
- # First try PyFITS for sanity
- testfile = os.path.join(datadir, "test"+tchar[i]+".fits")
- test_array = pyfits.getdata(testfile)
- np.testing.assert_array_equal(fits_ref_array.astype(types[i]), test_array,
+ array_type = types[i]
+
+ #
+ # Test input from a single external FITS image
+ #
+
+ # Read the reference image to from an externally-generated fits file
+ test_file = os.path.join(datadir, "test"+tchar[i]+".fits")
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_array,
err_msg="PyFITS failing to read reference image.")
- # Then use the Image methods... Note this also relies on the array look working too
- image_init_func = eval("galsim.Image"+tchar[i]) # Use handy eval() mimics use of ImageSIFD
- # First give ImageSIFD.read() a PyFITS PrimaryHDU
- hdu = pyfits.open(testfile)
- # NB: The returned test_image is not necessarily of type image_init_func!
- test_image = image_init_func.read(hdu)
- np.testing.assert_array_equal(fits_ref_array.astype(types[i]), test_image.array,
- err_msg="Image"+tchar[i]+".read() failed reading from PyFITS PrimaryHDU input.")
- # Then try an ImageSIFD.read() with the filename itself as input
- test_image = image_init_func.read(testfile)
- np.testing.assert_array_equal(fits_ref_array.astype(types[i]), test_image.array,
+
+ # Then use galsim fits.read function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_file)
+ test_image = galsim.fits.read(hdu)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_image.array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image = galsim.fits.read(test_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_image.array,
err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
- # TODO: test reading from an HDU list (e.g. for multi-extension FITS).
+
+ #
+ # Test full I/O on a single internally-generated FITS image
+ #
+
+ # Write the reference image to a fits file
+ ref_image = galsim.ImageView[array_type](ref_array.astype(array_type))
+ test_file = os.path.join(datadir, "test"+tchar[i]+"_internal.fits")
+ ref_image.write(test_file)
+
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_array,
+ err_msg="PyFITS failing to read reference image.")
+
+ # Then use galsim fits.read function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_file)
+ test_image = galsim.fits.read(hdu)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_image.array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image = galsim.fits.read(test_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_image.array,
+ err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
+
+ #
+ # Test input from an external multi-extension fits file
+ #
+
+ test_multi_file = os.path.join(datadir, "test_multi"+tchar[i]+".fits")
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_multi_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_array,
+ err_msg="PyFITS failing to read multi file.")
+
+ # Then use galsim fits.readMulti function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_multi_file)
+ test_image_list = galsim.fits.readMulti(hdu)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image_list = galsim.fits.readMulti(test_multi_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
+
+ #
+ # Test full I/O for an internally-generated multi-extension fits file
+ #
+
+ # Build a list of images with different values
+ image_list = []
+ for k in range(nimages):
+ image_list.append(ref_image + k)
+
+ # Write the list to a multi-extension fits file
+ test_multi_file = os.path.join(datadir, "test_multi"+tchar[i]+"_internal.fits")
+ galsim.fits.writeMulti(image_list,test_multi_file)
+
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_multi_file)
+ np.testing.assert_array_equal(ref_array.astype(types[i]), test_array,
+ err_msg="PyFITS failing to read multi file.")
+
+ # Then use galsim fits.readMulti function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_multi_file)
+ test_image_list = galsim.fits.readMulti(hdu)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image_list = galsim.fits.readMulti(test_multi_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
+
+ #
+ # Test input from an external fits data cube
+ #
+ test_cube_file = os.path.join(datadir, "test_cube"+tchar[i]+".fits")
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_cube_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]), test_array[k,:,:],
+ err_msg="PyFITS failing to read cube file.")
+
+ # Then use galsim fits.readCube function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_cube_file)
+ test_image_list = galsim.fits.readCube(hdu)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image_list = galsim.fits.readCube(test_cube_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
+
+ #
+ # Test full I/O for an internally-generated fits data cube
+ #
+
+ # Write the list to a fits data cube
+ test_cube_file = os.path.join(datadir, "test_cube"+tchar[i]+"_internal.fits")
+ galsim.fits.writeCube(image_list,test_cube_file)
+
+ # Check pyfits read for sanity
+ test_array = pyfits.getdata(test_cube_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]), test_array[k,:,:],
+ err_msg="PyFITS failing to read cube file.")
+
+ # Then use galsim fits.readCube function
+ # First version: use pyfits HDUList
+ hdu = pyfits.open(test_cube_file)
+ test_image_list = galsim.fits.readCube(hdu)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Failed reading from PyFITS PrimaryHDU input.")
+
+ # Second version: use file name
+ test_image_list = galsim.fits.readCube(test_cube_file)
+ for k in range(nimages):
+ np.testing.assert_array_equal((ref_array+k).astype(types[i]),
+ test_image_list[k].array,
+ err_msg="Image"+tchar[i]+".read() failed reading from string filename input.")
+
+
def test_Image_array_view():
"""Test that all four types of supported Images correctly provide a view on an input array.
Something went wrong with that request. Please try again.