Skip to content

Commit 69cab02

Browse files
briantobybriantoby
committed
first working HDF5-MAXIV importer
1 parent 9db660f commit 69cab02

File tree

1 file changed

+89
-89
lines changed

1 file changed

+89
-89
lines changed

‎GSASII/imports/G2pwd_HDF5.py‎

Lines changed: 89 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ class HDF5_Reader(G2obj.ImportPowderData):
4444
Any parameters placed in that file will override values set in the HDF5
4545
file.
4646
'''
47-
mode = None
47+
#mode = None
4848
def __init__(self):
4949
if h5py is None:
5050
self.UseReader = False
@@ -58,9 +58,12 @@ def __init__(self):
5858

5959
def ShowH5Element(self,obj,keylist):
6060
'''Format the contents of an HDF5 entry as a single line. Not used for
61-
reading files, only for use in :meth:`HDF5list`
61+
reading files, only used in :meth:`HDF5list`
6262
'''
6363
k = '/'.join(keylist)
64+
l = obj.get(k, getlink=True)
65+
if isinstance(l, h5py.ExternalLink):
66+
return f'link to file {l.filename}'
6467
try:
6568
typ = str(type(obj[k]))
6669
except:
@@ -75,15 +78,16 @@ def ShowH5Element(self,obj,keylist):
7578
return f'value={bool(obj[k][()])}'
7679
elif datfmt in ('<f8', 'uint8', 'int64', '<f4'): # scalar value or array of values
7780
try:
78-
return f'length {len(obj[k][()])}'
81+
len(obj[k][()])
82+
return f'array {obj[k].shape}'
7983
except:
8084
return f'value={obj[k][()]}'
8185
else:
8286
return f'dataset of type {repr(datfmt)}'
8387
elif ".Group'" in typ:
8488
return "(group)"
8589
else:
86-
return f'{prefix}/{i} is {type(obj[k])}'
90+
return f'type is {type(obj[k])}'
8791

8892
def RecurseH5Element(self,obj,prefix=[]):
8993
'''Returns a list of entries of all keys in the HDF5 file
@@ -94,7 +98,7 @@ def RecurseH5Element(self,obj,prefix=[]):
9498
* entry 0 is the top-level keys (/a, /b,...),
9599
* entry 1 has the first level keys (/a/c /a/d, /b/d, /b/e,...)
96100
* ...
97-
Not used for reading files, only for use in :meth:`HDF5list`
101+
Not used for reading files, used only in :meth:`HDF5list`
98102
'''
99103
try:
100104
self.HDF5entries
@@ -106,13 +110,16 @@ def RecurseH5Element(self,obj,prefix=[]):
106110
for i in obj:
107111
nextprefix = prefix+[i]
108112
self.HDF5entries[depth].append(nextprefix)
113+
# check for link objects
114+
l = obj.get(i, getlink=True)
115+
if isinstance(l, h5py.ExternalLink): continue
109116
try:
110117
typ = str(type(obj[i]))
111118
except:
112119
print(f'**Error** with key {prefix}/{i}')
113120
continue
114121
if ".Group'" in typ:
115-
t = f'{prefix}/{i}'
122+
#t = f'{prefix}/{i}'
116123
#print(f'\n{nextprefix} contents {(60-len(t))*'='}')
117124
self.RecurseH5Element(obj[i],nextprefix)
118125
return self.HDF5entries
@@ -137,10 +144,12 @@ def HDF5list(self, filename):
137144
m = max(m,len('/'.join(k)))
138145
for k in j:
139146
lbl = self.ShowH5Element(fp,k)
147+
if '\n' in lbl:
148+
lbl = '; '.join(lbl.split('\n'))
140149
if len(lbl) > 50:
141150
lbl = lbl[:50] + '...'
142-
if '\n' in lbl:
143-
lbl = lbl.split()[0] + '...'
151+
# if '\n' in lbl:
152+
# lbl = lbl.split()[0] + '...'
144153
if lbl != '(group)': strings.append(f"{'/'.join(k):{m}s} {lbl}")
145154
with open(filename+'_contents.txt', 'w') as fp:
146155
for i in strings: fp.write(f'{i}\n')
@@ -149,17 +158,19 @@ def ContentsValidator(self, filename):
149158
'''Test if valid by seeing if the HDF5 library recognizes the file.
150159
Then get file type (currently MAX IV NeXus/NXazint[12]d only)
151160
'''
152-
from .. import GSASIIpath
161+
#from .. import GSASIIpath
153162
try:
154163
fp = h5py.File(filename, 'r')
155164
if 'entry' in fp: # NeXus
165+
#self.HDF5entries = []
166+
#self.HDF5list(filename)
156167
if 'definition' in fp['/entry']:
157168
# MAX IV NXazint1d file
158169
if fp['/entry/definition'][()].decode() == 'NXazint1d':
159170
return True
160171
# MAX IV NXazint1d file
161-
if fp['/entry/definition'][()].decode() == 'NXazint2d':
162-
return True
172+
#if fp['/entry/definition'][()].decode() == 'NXazint2d':
173+
# return True
163174
except IOError:
164175
return False
165176
finally:
@@ -182,9 +193,6 @@ def Reader(self, filename, ParentFrame=None, **kwarg):
182193
self.blknum = 0
183194
else:
184195
self.blknum = min(self.selections)
185-
try:
186-
self.mode = None
187-
fp = h5py.File(filename, 'r')
188196
try:
189197
fp = h5py.File(filename, 'r')
190198
if 'entry' in fp: # NeXus
@@ -211,103 +219,95 @@ def readNXazint1d(self, filename, fpbuffer={}):
211219
'''Read HDF5 file in NeXus as produced by MAX IV as a NXazint1d
212220
see https://nxazint-hdf5-nexus-3229ecbd09ba8a773fbbd8beb72cace6216dfd5063e1.gitlab-pages.esrf.fr/classes/contributed_definitions/NXazint1d.html
213221
'''
214-
self.instmsg = 'HDF file'
222+
#self.instmsg = 'HDF file'
215223
doread = False # has the file already been read into a buffer?
216-
for i in ('blkmap','intenArr_blknum')+self.midassections:
224+
arrays = ('entry/data/radial_axis','entry/data/I')
225+
floats = ('entry/instrument/monochromator/wavelength',
226+
'entry/reduction/input/polarization_factor')
227+
strings = ('entry/instrument/source/name','entry/reduction/input/unit')
228+
for i in arrays+floats+strings:
217229
if i not in fpbuffer:
218230
doread = True
219231
break
220-
else: # do we have the right section buffered?
221-
doread = fpbuffer['intenArr_blknum'] != self.blknum
222-
223232
if doread: # read into buffer
224233
try:
225234
fp = h5py.File(filename, 'r')
226-
if 'blkmap' not in fpbuffer:
227-
fpbuffer['blkmap'] = list(fp.get('OmegaSumFrame').keys())
228-
if 'REtaMap' not in fpbuffer:
229-
fpbuffer['REtaMap'] = np.array(fp.get('REtaMap'))
230-
if 'intenArr' not in fpbuffer or fpbuffer.get('intenArr_blknum',-1) != self.blknum:
231-
fpbuffer['intenArr'] = np.array(fp.get('OmegaSumFrame').get(
232-
fpbuffer['blkmap'][self.blknum]))
233-
fpbuffer['intenArr_blknum'] = self.blknum
234-
self.azmcnt = -1
235-
if 'Omegas' not in fpbuffer:
236-
fpbuffer['Omegas'] = np.array(fp.get('Omegas'))
235+
for i in arrays:
236+
fpbuffer[i] = np.array(fp.get(i))
237+
for i in floats:
238+
fpbuffer[i] = float(fp[i][()])
239+
for i in strings:
240+
fpbuffer[i] = fp[i][()].decode()
241+
if fpbuffer['entry/reduction/input/unit'] != '2th':
242+
print('NXazint1d HDF5 file has units',fpbuffer['entry/reduction/input/unit'])
243+
self.errors = 'NXazint1d only can be read with 2th units'
244+
return False
245+
if self.selections is None or len(self.selections) == 0:
246+
self.blknum = 0
247+
else:
248+
self.blknum = min(self.selections)
237249
except IOError:
238250
print ('cannot open file '+ filename)
239251
return False
240252
finally:
241253
fp.close()
242-
# get overriding sample & instrument parameters
243-
fpbuffer['sampleprm'] = {}
244-
samplefile = os.path.splitext(filename)[0] + '.samprm'
245-
if os.path.exists(samplefile):
246-
fp = open(samplefile,'r')
247-
S = fp.readline()
248-
while S:
249-
if not S.strip().startswith('#'):
250-
[item,val] = S[:-1].split(':')
251-
fpbuffer['sampleprm'][item.strip("'")] = eval(val)
252-
S = fp.readline()
253-
fp.close()
254-
fpbuffer['instprm'] = {}
255-
instfile = os.path.splitext(filename)[0] + '.instprm'
256-
if os.path.exists(instfile):
257-
self.instmsg = 'HDF and .instprm files'
258-
fp = open(instfile,'r')
259-
S = fp.readline()
260-
while S:
261-
if not S.strip().startswith('#'):
262-
[item,val] = S[:-1].split(':')
263-
fpbuffer['instprm'][item.strip("'")] = eval(val)
264-
S = fp.readline()
265-
fp.close()
266-
# look for a non-empty scan (lineout)
267-
use = [0]
268-
while sum(use) == 0 and self.azmcnt < fpbuffer['intenArr'].shape[1]:
269-
self.azmcnt += 1
270-
if self.azmcnt >= fpbuffer['intenArr'].shape[1]:
271-
return False
272-
use = fpbuffer['REtaMap'][3,:,self.azmcnt] != 0
273-
254+
self.numbanks=len(fpbuffer['entry/data/I'])
255+
# # get overriding sample & instrument parameters
256+
# fpbuffer['sampleprm'] = {}
257+
# samplefile = os.path.splitext(filename)[0] + '.samprm'
258+
# if os.path.exists(samplefile):
259+
# fp = open(samplefile,'r')
260+
# S = fp.readline()
261+
# while S:
262+
# if not S.strip().startswith('#'):
263+
# [item,val] = S[:-1].split(':')
264+
# fpbuffer['sampleprm'][item.strip("'")] = eval(val)
265+
# S = fp.readline()
266+
# fp.close()
267+
# fpbuffer['instprm'] = {}
268+
# instfile = os.path.splitext(filename)[0] + '.instprm'
269+
# if os.path.exists(instfile):
270+
# self.instmsg = 'HDF and .instprm files'
271+
# fp = open(instfile,'r')
272+
# S = fp.readline()
273+
# while S:
274+
# if not S.strip().startswith('#'):
275+
# [item,val] = S[:-1].split(':')
276+
# fpbuffer['instprm'][item.strip("'")] = eval(val)
277+
# S = fp.readline()
278+
# fp.close()
274279
# now transfer information into current histogram
275-
self.pwdparms['Instrument Parameters'] = [
276-
{'Type': ['PXC', 'PXC', False]},
277-
{}]
278-
inst = {}
279-
inst.update(instprmList)
280-
inst.update(fpbuffer['instprm'])
281-
for key,val in inst.items():
282-
self.pwdparms['Instrument Parameters'][0][key] = [val,val,False]
283-
samp = {}
284-
samp.update(sampleprmList)
285-
samp.update(fpbuffer['sampleprm'])
286-
for key,val in samp.items():
287-
self.Sample[key] = val
288-
self.numbanks=len(fpbuffer['blkmap'])
289-
x = fpbuffer['REtaMap'][1,:,self.azmcnt][use]
290-
y = fpbuffer['intenArr'][:,self.azmcnt][use]
291-
w = np.nan_to_num(1/y) # this is probably not correct
292-
eta = np.average(fpbuffer['REtaMap'][2,:,self.azmcnt][use])
293-
self.pwdparms['Instrument Parameters'][0]['Azimuth'] = [90-eta,90-eta,False]
294-
self.pwdparms['Instrument Parameters'][0]['Bank'] = [self.azmcnt,self.azmcnt,False]
280+
#self.pwdparms['Instrument Parameters'] = [
281+
# {'Type': ['PXC', 'PXC', False]},
282+
# {}]
283+
# inst = {}
284+
# inst.update(instprmList)
285+
# inst.update(fpbuffer['instprm'])
286+
# for key,val in inst.items():
287+
# self.pwdparms['Instrument Parameters'][0][key] = [val,val,False]
288+
# samp = {}
289+
# samp.update(sampleprmList)
290+
# samp.update(fpbuffer['sampleprm'])
291+
# for key,val in samp.items():
292+
# self.Sample[key] = val
293+
x = fpbuffer['entry/data/radial_axis']
294+
y = fpbuffer['entry/data/I'][self.blknum]
295+
w = np.nan_to_num(1/y) # this is not correct
296+
#self.pwdparms['Instrument Parameters'][0]['Azimuth'] = [90-eta,90-eta,False]
297+
#self.pwdparms['Instrument Parameters'][0]['Bank'] = [self.blknum,self.blknum,False]
295298
# self.Sample['Gonio. radius'] = float(S.split('=')[1])
296299
# self.Sample['Omega'] = float(S.split('=')[1])
297300
# self.Sample['Chi'] = float(S.split('=')[1])
298-
self.Sample['Phi'] = Omega = fpbuffer['Omegas'][self.blknum]
301+
#self.Sample['Phi'] = Omega = fpbuffer['Omegas'][self.blknum]
299302
self.powderdata = [x,y,w,np.zeros_like(x),np.zeros_like(x),np.zeros_like(x)]
300303
#self.comments = comments[selblk]
301304
self.powderentry[0] = filename
302305
#self.powderentry[1] = Pos # position offset (never used, I hope)
303306
self.powderentry[2] = self.blknum # bank number
304-
self.idstring = f'{os.path.split(filename)[1][:10]} omega={Omega} eta={eta}'
305-
# if GSASIIpath.GetConfigValue('debug'): print(
306-
# f'Read entry #{self.azmcnt} img# {self.blknum} from file {filename}')
307-
# are there more lineouts after this one in current image to read?
308-
self.repeat = sum(sum(fpbuffer['REtaMap'][3,:,self.azmcnt+1:])) != 0
309-
if self.repeat: return True
307+
self.idstring = f'#{self.blknum} {os.path.split(filename)[1][:60]}'
308+
self.instdict['wave'] = fpbuffer['entry/instrument/monochromator/wavelength']
310309
# if not, are there more [selected] images that after this to be read?
310+
self.repeat = False
311311
if self.blknum < self.numbanks-1:
312312
if self.selections is None or len(self.selections) == 0:
313313
self.blknum += 1

0 commit comments

Comments
 (0)