-
Notifications
You must be signed in to change notification settings - Fork 17
/
lp.py
212 lines (196 loc) · 6.78 KB
/
lp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
"""
Set of tools to read and write 'La Palma' cubes
"""
import os
import numpy as np
def make_header(image):
''' Creates header for La Palma images. '''
from struct import pack
ss = image.shape
# only 2D or 3D arrays
if len(ss) not in [2, 3]:
raise IndexError('make_header: input array must be 2D or 3D, got %iD'
% len(ss))
dtypes = {'int8': ['(byte)', 1], 'int16': ['(integer)', 2],
'int32': ['(long)', 3], 'float32': ['(float)', 4]}
if str(image.dtype) not in dtypes:
raise ValueError('make_header: array type' +
' %s not supported, must be one of %s' %
(image.dtype, list(dtypes.keys())))
sdt = dtypes[str(image.dtype)]
header = ' datatype=%s %s, dims=%i, nx=%i, ny=%i' % \
(sdt[1], sdt[0], len(ss), ss[0], ss[1])
if len(ss) == 3:
header += ', nt=%i' % (ss[2])
# endianess
if pack('@h', 1) == pack('<h', 1):
header += ', endian=l'
else:
header += ', endian=b'
return header
def writeto(filename, image, extraheader='', dtype=None, verbose=False,
append=False):
'''Writes image into cube, La Palma format. Analogous to IDL's lp_write.'''
# Tiago notes: seems to have problems with 2D images, but not sure if that
# even works in IDL's routines...
if not os.path.isfile(filename):
append = False
# use dtype from array, if none is specified
if dtype is None:
dtype = image.dtype
image = image.astype(dtype)
if append:
# check if image sizes/types are consistent with file
sin, t, h = getheader(filename)
if sin[:2] != image.shape[:2]:
raise IOError('writeto: trying to write' +
' %s images, but %s has %s images!' %
(repr(image.shape[:2]), filename, repr(sin[:2])))
if np.dtype(t) != image.dtype:
raise IOError('writeto: trying to write' +
' %s type images, but %s nas %s images' %
(image.dtype, filename, np.dtype(t)))
# add the nt of current image to the header
hloc = h.lower().find('nt=')
new_nt = str(sin[-1] + image.shape[-1])
header = h[:hloc + 3] + new_nt + h[hloc + 3 + len(str(sin[-1])):]
else:
header = make_header(image)
if extraheader:
header += ' : ' + extraheader
# convert string to [unsigned] byte array
hh = np.zeros(512, dtype='uint8')
for i, ss in enumerate(header):
hh[i] = ord(ss)
# write header to file
file_arr = np.memmap(filename, dtype='uint8', mode=append and 'r+' or 'w+',
shape=(512,))
file_arr[:512] = hh[:]
del file_arr
# offset if appending
apoff = append and np.prod(sin) * image.dtype.itemsize or 0
# write array to file
file_arr = np.memmap(filename, dtype=dtype, mode='r+', order='F',
offset=512 + apoff, shape=image.shape)
file_arr[:] = image[:]
del file_arr
if verbose:
if append:
print(('Appended %s %s array into %s.' % (image.shape, dtype,
filename)))
else:
print(('Wrote %s, %s array of shape %s' % (filename, dtype,
image.shape)))
return
def writeheader(filename, header):
"""
Writes header (proper format, from make_header) into existing filename.
"""
# convert string to [unsigned] byte array
hh = np.zeros(512, dtype='uint8')
for i, ss in enumerate(header):
hh[i] = ord(ss)
# write header to file
file_arr = np.memmap(filename, dtype='uint8', mode='r+', shape=(512,))
file_arr[:512] = hh[:]
del file_arr
return
def getheader(filename):
"""
Reads header from La Palma format cube.
Parameters
----------
filename : str
File to read header from.
Returns
-------
result : list
List with shape tuple (nx, ny [, nt]),
datatype (with endianness), header string.
"""
# read header and convert to string
h = np.fromfile(filename, dtype='uint8', count=512)
header = ''
for s in h[h > 0]:
header += chr(s)
# start reading at 'datatype'
hd = header[header.lower().find('datatype'):]
hd = hd.split(':')[0].replace(',', ' ').split()
# Types: uint8 int16 int32 float32
typelist = ['u1', 'i2', 'i4', 'f4']
# extract datatype
try:
dtype = typelist[int(hd[0].split('=')[1]) - 1]
except:
print(header)
raise IOError('getheader: datatype invalid or missing')
# extract endianness
try:
if hd[-1].split('=')[0].lower() != 'endian':
raise IndexError()
endian = hd[-1].split('=')[1]
except IndexError:
print(header)
raise IOError('getheader: endianess missing.')
if endian.lower() == 'l':
dtype = '<' + dtype
else:
dtype = '>' + dtype
# extract dims
try:
if hd[2].split('=')[0].lower() != 'dims':
raise IndexError()
dims = int(hd[2].split('=')[1])
if dims not in [2, 3]:
raise ValueError('Invalid dims=%i (must be 2 or 3)' % dims)
except IndexError:
print(header)
raise IOError('getheader: dims invalid or missing.')
try:
if hd[3].split('=')[0].lower() != 'nx':
raise IndexError()
nx = int(hd[3].split('=')[1])
except:
print(header)
raise IOError('getheader: nx invalid or missing.')
try:
if hd[4].split('=')[0].lower() != 'ny':
raise IndexError()
ny = int(hd[4].split('=')[1])
except:
print(header)
raise IOError('getheader: ny invalid or missing.')
if dims == 3:
try:
if hd[5].split('=')[0].lower() != 'nt':
raise IndexError()
nt = int(hd[5].split('=')[1])
except:
print(header)
raise IOError('getheader: nt invalid or missing.')
shape = (nx, ny, nt)
else:
shape = (nx, ny)
return [shape, dtype, header]
def getdata(filename, rw=False, verbose=False):
"""
Reads La Palma format cube into a memmap object.
Parameters
----------
filename : str
File to read.
rw : bool, optional
If True, then any change to the data will be written to file.
verbose : bool, optional
If True, will print out additional information.
Returns
-------
data : ndarray
Numpy memmap object with data.
"""
sh, dt, header = getheader(filename)
if verbose:
print(('Reading %s...\n%s' % (filename, header)))
mode = ['c', 'r+']
return np.memmap(filename, mode=mode[rw], shape=sh, dtype=dt, order='F',
offset=512)