Merge pull request #617 from astrofrog/big-endian

Add support for memmap=False and big endian CASA images

spectral_cube/io/casa_dask.py

@@ -114,6 +114,9 @@ def casa_image_dask_reader(imagename, memmap=True, mask=False):
     # tb.close()
     dminfo = getdminfo(str(imagename))
 
+    # Determine whether file is big endian
+    big_endian = dminfo['*1']['BIGENDIAN']
+
     # chunkshape defines how the chunks (array subsets) are written to disk
     chunkshape = tuple(dminfo['*1']['SPEC']['DEFAULTTILESHAPE'])
     chunksize = np.product(chunkshape)
@@ -141,10 +144,16 @@ def casa_image_dask_reader(imagename, memmap=True, mask=False):
                              "expected {1}".format(filesize, expected))
     else:
         if filesize == totalsize * 4:
-            dtype = np.float32
+            if big_endian:
+                dtype = '>f4'
+            else:
+                dtype = '<f4'
             itemsize = 4
         elif filesize == totalsize * 8:
-            dtype = np.float64
+            if big_endian:
+                dtype = '>f8'
+            else:
+                dtype = '<f8'
             itemsize = 8
         else:
             raise ValueError("Unexpected file size for data, found {0} but "
@@ -165,11 +174,11 @@ def casa_image_dask_reader(imagename, memmap=True, mask=False):
             full_array = np.unpackbits(full_array, bitorder='little').astype(np.bool_)
             ceil_chunksize = int(ceil(chunksize / 8)) * 8
             chunks = [full_array[ii*ceil_chunksize:(ii+1)*ceil_chunksize][:chunksize].reshape(chunkshape, order='F').T
-                    for ii in range(nchunks)]
+                      for ii in range(nchunks)]
         else:
             full_array = np.fromfile(img_fn, dtype=dtype)
             chunks = [full_array[ii*chunksize:(ii+1)*chunksize].reshape(chunkshape, order='F').T
-                    for ii in range(nchunks)]
+                      for ii in range(nchunks)]
 
     # convert all chunks to dask arrays - and set name and meta appropriately
     # to prevent dask trying to access the data to determine these

spectral_cube/io/casa_image.py

@@ -38,7 +38,7 @@ def is_casa_image(origin, filepath, fileobj, *args, **kwargs):
     return filepath is not None and filepath.lower().endswith('.image')
 
 
-def load_casa_image(filename, skipdata=False,
+def load_casa_image(filename, skipdata=False, memmap=True,
                     skipvalid=False, skipcs=False, target_cls=None, **kwargs):
     """
     Load a cube (into memory?) from a CASA image. By default it will transpose
@@ -53,14 +53,14 @@ def load_casa_image(filename, skipdata=False,
 
     # read in the data
     if not skipdata:
-        data = casa_image_dask_reader(filename)
+        data = casa_image_dask_reader(filename, memmap=memmap)
 
     # CASA stores validity of data as a mask
     if skipvalid:
         valid = None
     else:
         try:
-            valid = casa_image_dask_reader(filename, mask=True)
+            valid = casa_image_dask_reader(filename, memmap=memmap, mask=True)
         except FileNotFoundError:
             valid = None
 

spectral_cube/io/casa_low_level_io.py

@@ -224,9 +224,11 @@ def read_table(f, image_path):
         raise NotImplementedError('Support for {0} version {1} not implemented'.format(stype, sversion))
 
     nrow = read_int32(f)
-    unknown = read_int32(f)  # noqa
+    fmt = read_int32(f)  # noqa
     name = read_string(f)  # noqa
 
+    big_endian = fmt == 0  # noqa
+
     table_desc = read_table_desc(f, nrow, image_path)
 
     return table_desc
@@ -317,7 +319,7 @@ def read_tiled_st_man(f):
     st_man = {}
     st_man['SPEC'] = {}
 
-    big_endian = f.read(1)  # noqa
+    st_man['BIGENDIAN'] = f.read(1) == b'\x01'  # noqa
 
     seqnr = read_int32(f)
     if seqnr != 0:

spectral_cube/io/tests/test_casa_low_level_io.py

@@ -55,6 +55,9 @@ def test_getdminfo(tmp_path, shape):
 
     actual = getdminfo(filename)
 
+    # We include information about endian-ness in the dminfo but CASA doesn't
+    actual['*1'].pop('BIGENDIAN')
+
     # The easiest way to compare the output is simply to compare the output
     # from pformat (checking for dictionary equality doesn't work because of
     # the Numpy arrays inside).

spectral_cube/tests/data/basic_bigendian.image/logtable/table.info

@@ -0,0 +1,4 @@
+Type = Log message
+SubType = 
+
+Repository for software-generated logging messages

spectral_cube/tests/data/basic_bigendian.image/mask0/table.f0_TSM0

@@ -0,0 +1 @@
+���������������

spectral_cube/tests/data/basic_bigendian.image/mask0/table.info

@@ -0,0 +1,3 @@
+Type = Paged Array
+SubType = 
+

spectral_cube/tests/data/basic_bigendian.image/table.info

@@ -0,0 +1,3 @@
+Type = Image
+SubType = 
+

spectral_cube/tests/test_casafuncs.py

@@ -2,6 +2,7 @@
 
 import os
 import shutil
+from itertools import product
 
 import pytest
 import numpy as np
@@ -72,14 +73,19 @@ def filename(request):
     return request.getfixturevalue(request.param)
 
 
-def test_casa_read_basic():
+@pytest.mark.parametrize(('memmap', 'bigendian'), product((False, True), (False, True)))
+def test_casa_read_basic(memmap, bigendian):
 
     # Check that SpectralCube.read works for an example CASA dataset stored
     # in the tests directory. This test should NOT require CASA, whereas a
     # number of tests below require CASA to generate test datasets. The present
     # test is to ensure CASA is not required for reading.
 
-    cube = SpectralCube.read(os.path.join(DATA, 'basic.image'))
+    if bigendian:
+        cube = SpectralCube.read(os.path.join(DATA, 'basic_bigendian.image'), memmap=memmap)
+    else:
+        cube = SpectralCube.read(os.path.join(DATA, 'basic.image'), memmap=memmap)
+
     assert cube.shape == (3, 4, 5)
     assert_allclose(cube.wcs.pixel_to_world_values(1, 2, 3),
                     [2.406271e+01, 2.993521e+01, 1.421911e+09])