source: trunk/GSASIIfiles.py

Last change on this file was 5112, checked in by vondreele, 10 months ago

implement energy dispersive x-ray data as a new type 'PXE'; assume Gaussian peak shape only; no sample broadening considered. Peak resolution is only ~0.5% not good enough to see sample broadening
import, plotting, peak fitting, indexing & cell refinement from peak positions "complete"
data is old (from ~20 yrs ago) so no idea about modern ED data.
Start on use in Pawley/LeBail? refinement - Rietveld excluded; currently not working
Some work on ISODISTORT implementation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 45.6 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-12-15 22:56:41 +0000 (Wed, 15 Dec 2021) $
4# $Author: vondreele $
5# $Revision: 5112 $
6# $URL: trunk/GSASIIfiles.py $
7# $Id: GSASIIfiles.py 5112 2021-12-15 22:56:41Z vondreele $
8########### SVN repository information ###################
9'''
10*GSASIIfiles: data (non-GUI) I/O routines*
11==========================================
12
13Module with miscellaneous routines for input and output from files.
14
15This module should not contain any references to wxPython so that it
16can be imported for scriptable use or potentially on clients where
17wx is not installed.
18
19Future refactoring: This module and GSASIIIO.py needs some work to
20move non-wx routines here. It may will likely make sense to rename the module(s)
21at that point.
22'''
23from __future__ import division, print_function
24import platform
25import os
26import sys
27import glob
28import inspect
29import numpy as np
30
31import GSASIIpath
32GSASIIpath.SetVersionNumber("$Revision: 5112 $")
33
34# N.B. This is duplicated in G2IO
35def sfloat(S):
36    'Convert a string to float. An empty field or a unconvertable value is treated as zero'
37    if S.strip():
38        try:
39            return float(S)
40        except ValueError:
41            pass
42    return 0.0
43
44G2printLevel = 'all'
45'''This defines the level of output from calls to :func:`G2Print`,
46which should  be used in place of print() within this module.
47Settings for this are 'all', 'warn', 'error' or 'none'. Also see:
48:func:`G2Print` and :func:`G2SetPrintLevel`.
49'''
50
51def G2SetPrintLevel(level):
52    '''Set the level of output from calls to :func:`G2Print`, which should
53    be used in place of print() within GSASII. Settings for the mode are
54    'all', 'warn', 'error' or 'none'
55   
56    :param str level: a string used to set the print level, which may be
57      'all', 'warn', 'error' or 'none'.
58      Note that capitalization and extra letters in level are ignored, so
59      'Warn', 'warnings', etc. will all set the mode to 'warn'
60    '''
61    global G2printLevel
62    for mode in  'all', 'warn', 'error', 'none':
63        if mode in level.lower():
64            G2printLevel = mode
65            return
66    else:
67        G2Print('G2SetPrintLevel Error: level={} cannot be interpreted.',
68                    'Use all, warn, error or none.')
69       
70def find(name, path):
71    '''find 1st occurance of file in path
72    '''
73    for root, dirs, files in os.walk(path):
74        if name in files:
75            return os.path.join(root, name)
76       
77def G2Print(*args,**kwargs):
78    '''Print with filtering based level of output (see :func:`G2SetPrintLevel`).
79    Use G2Print() as replacement for print().
80
81    :param str mode: if specified, this should contain the mode for printing
82      ('error', 'warn' or anything else). If not specified, the first argument
83      of the print command (args[0]) should contain the string 'error' for
84      error messages and 'warn' for warning messages
85      (capitalization and additional letters ignored.)
86    '''
87    if G2printLevel == 'none': return
88    if kwargs.get('mode') is None:
89        testStr = str(args[0]).lower()
90    else:
91        testStr = kwargs['mode'][:].lower()
92        del kwargs['mode'] 
93    level = 2
94    for i,mode in enumerate(('error', 'warn')):
95        if mode in testStr:
96            level = i
97            break
98    if G2printLevel == 'error' and level > 0: return
99    if G2printLevel == 'warn' and level > 1: return
100    print(*args,**kwargs)
101
102def get_python_versions(packagelist):
103    versions = [['Python', sys.version.split()[0]]]
104    for pack in packagelist:
105        try:
106            versions.append([pack.__name__, pack.__version__])
107        except:
108            pass
109    versions.append(['Platform',
110                     sys.platform + ' ' + platform.architecture()[0] +
111                     ' ' + platform.machine()])
112    return versions
113
114def makeInstDict(names,data,codes):
115    inst = dict(zip(names,zip(data,data,codes)))
116    for item in inst:
117        inst[item] = list(inst[item])
118    return inst
119
120def SetPowderInstParms(Iparm, rd):
121    '''extracts values from instrument parameters in rd.instdict
122    or in array Iparm.
123    Create and return the contents of the instrument parameter tree entry.
124    '''
125    Irads = {0:' ',1:'CrKa',2:'FeKa',3:'CuKa',4:'MoKa',5:'AgKa',6:'TiKa',7:'CoKa'}
126    DataType = Iparm['INS   HTYPE '].strip()[:3]  # take 1st 3 chars
127    # override inst values with values read from data file
128    Bank = rd.powderentry[2]    #should be used in multibank iparm files
129    if rd.instdict.get('type'):
130        DataType = rd.instdict.get('type')
131    data = [DataType,]
132    instname = Iparm.get('INS  1INAME ')
133    irad = int(Iparm.get('INS  1 IRAD ','0'))
134    if instname:
135        rd.Sample['InstrName'] = instname.strip()
136    if 'C' in DataType:
137        wave1 = None
138        wave2 = 0.0
139        if rd.instdict.get('wave'):
140            wl = rd.instdict.get('wave')
141            wave1 = wl[0]
142            if len(wl) > 1: wave2 = wl[1]
143        s = Iparm['INS  1 ICONS']
144        if not wave1:
145            wave1 = sfloat(s[:10])
146            wave2 = sfloat(s[10:20])
147        v = (wave1,wave2,
148             sfloat(s[20:30])/100.,sfloat(s[55:65]),sfloat(s[40:50])) #get lam1, lam2, zero, pola & ratio
149        if not v[1]:
150            names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth']
151            v = (v[0],v[2],v[4])
152            codes = [0,0,0,0,0]
153            rd.Sample.update({'Type':'Debye-Scherrer','Absorption':[0.,False],'DisplaceX':[0.,False],'DisplaceY':[0.,False]})
154        else:
155            names = ['Type','Lam1','Lam2','Zero','I(L2)/I(L1)','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth']
156            codes = [0,0,0,0,0,0,0]
157            rd.Sample.update({'Type':'Bragg-Brentano','Shift':[0.,False],'Transparency':[0.,False],
158                'SurfRoughA':[0.,False],'SurfRoughB':[0.,False]})
159        data.extend(v)
160        if 'INS  1PRCF  ' in Iparm:
161            v1 = Iparm['INS  1PRCF  '].split()
162            v = Iparm['INS  1PRCF 1'].split()
163            data.extend([float(v[0]),float(v[1]),float(v[2])])                  #get GU, GV & GW - always here
164            azm = float(Iparm.get('INS  1DETAZM','0.0'))
165            v = Iparm['INS  1PRCF 2'].split()
166            if v1[0] == 3:
167                data.extend([float(v[0]),float(v[1]),0.0,float(v[2])+float(v[3],azm)])  #get LX, LY, Z, S+H/L & azimuth
168            else:
169                data.extend([0.0,0.0,0.0,0.002,azm])                                      #OK defaults if fxn #3 not 1st in iprm file
170        else:
171            v1 = Iparm['INS  1PRCF1 '].split()
172            v = Iparm['INS  1PRCF11'].split()
173            data.extend([float(v[0]),float(v[1]),float(v[2])])                  #get GU, GV & GW - always here
174            azm = float(Iparm.get('INS  1DETAZM','0.0'))
175            v = Iparm['INS  1PRCF12'].split()
176            if v1[0] == 3:
177                data.extend([float(v[0]),float(v[1]),0.0,float(v[2])+float(v[3],azm)])  #get LX, LY, Z, S+H/L & azimuth
178            else:
179                data.extend([0.0,0.0,0.0,0.002,azm])                                      #OK defaults if fxn #3 not 1st in iprm file
180        codes.extend([0,0,0,0,0,0,0])
181        Iparm1 = makeInstDict(names,data,codes)
182        Iparm1['Source'] = [Irads[irad],Irads[irad]]
183        Iparm1['Bank'] = [Bank,Bank,0]
184        return [Iparm1,{}]
185    elif 'T' in DataType:
186        names = ['Type','fltPath','2-theta','difC','difA', 'difB','Zero','alpha','beta-0','beta-1',
187            'beta-q','sig-0','sig-1','sig-2','sig-q', 'X','Y','Z','Azimuth',]
188        codes = [0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,]
189        azm = 0.
190        if 'INS  1DETAZM' in Iparm:
191            azm = float(Iparm['INS  1DETAZM'])
192        rd.Sample['Azimuth'] = azm
193        fltPath0 = 20.                      #arbitrary
194        if 'INS   FPATH1' in Iparm:
195            s = Iparm['INS   FPATH1'].split()
196            fltPath0 = sfloat(s[0])
197        if 'INS  1BNKPAR' not in Iparm:     #bank missing from Iparm file
198            return []
199        s = Iparm['INS  1BNKPAR'].split()
200        fltPath1 = sfloat(s[0])
201        data.extend([fltPath0+fltPath1,])               #Flight path source-sample-detector
202        data.extend([sfloat(s[1]),])               #2-theta for bank
203        s = Iparm['INS  1 ICONS'].split()
204        data.extend([sfloat(s[0]),sfloat(s[1]),0.0,sfloat(s[2])])    #difC,difA,difB,Zero
205        if 'INS  1PRCF  ' in Iparm:
206            s = Iparm['INS  1PRCF  '].split()
207            pfType = int(s[0])
208            s = Iparm['INS  1PRCF 1'].split()
209            if abs(pfType) == 1:
210                data.extend([sfloat(s[1]),sfloat(s[2]),sfloat(s[3])]) #alpha, beta-0, beta-1
211                s = Iparm['INS  1PRCF 2'].split()
212                data.extend([0.0,0.0,sfloat(s[1]),sfloat(s[2]),0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
213            elif abs(pfType) in [3,4,5]:
214                data.extend([sfloat(s[0]),sfloat(s[1]),sfloat(s[2])]) #alpha, beta-0, beta-1
215                if abs(pfType) == 4:
216                    data.extend([0.0,0.0,sfloat(s[3]),0.0,0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
217                else:
218                    s = Iparm['INS  1PRCF 2'].split()
219                    data.extend([0.0,0.0,sfloat(s[0]),sfloat(s[1]),0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
220            elif abs(pfType) == 2:
221                print('''***WARNING gsas profile function #2 does not give valid GSAS-II diffractometer/profile coefficients ***
222                you should request GSAS-II instparm file from Instrument responsible''')
223                data.extend([sfloat(s[1]),0.0,1./sfloat(s[3])]) #alpha, beta-0, beta-1
224                data.extend([0.0,0.0,sfloat(s[1]),0.0,0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
225        else:
226            s = Iparm['INS  1PRCF1 '].split()
227            pfType = int(s[0])
228            s = Iparm['INS  1PRCF11'].split()
229            if abs(pfType) == 1:
230                data.extend([sfloat(s[1]),sfloat(s[2]),sfloat(s[3])]) #alpha, beta-0, beta-1
231                s = Iparm['INS  1PRCF12'].split()
232                data.extend([0.0,0.0,sfloat(s[1]),sfloat(s[2]),0.0,0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
233            elif abs(pfType) in [3,4,5]:
234                data.extend([sfloat(s[0]),sfloat(s[1]),sfloat(s[2])]) #alpha, beta-0, beta-1
235                if abs(pfType) == 4:
236                    data.extend([0.0,0.0,sfloat(s[3]),0.0,0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
237                else:
238                    s = Iparm['INS  1PRCF12'].split()
239                    data.extend([0.0,0.0,sfloat(s[0]),sfloat(s[1]),0.0,0.0,0.0,0.0,azm])    #beta-q, sig-0, sig-1, sig-2, sig-q, X, Y, Z
240        Inst1 = makeInstDict(names,data,codes)
241        Inst1['Bank'] = [Bank,Bank,0]
242        Inst2 = {}
243        if pfType < 0:
244            Ipab = 'INS  1PAB'+str(-pfType)
245            Npab = int(Iparm[Ipab+'  '].strip())
246            Inst2['Pdabc'] = []
247            for i in range(Npab):
248                k = Ipab+str(i+1).rjust(2)
249                s = Iparm[k].split()
250                Inst2['Pdabc'].append([float(t) for t in s])
251            Inst2['Pdabc'] = np.array(Inst2['Pdabc'])
252            Inst2['Pdabc'].T[3] += Inst2['Pdabc'].T[0]*Inst1['difC'][0] #turn 3rd col into TOF
253        if 'INS  1I ITYP' in Iparm:
254            s = Iparm['INS  1I ITYP'].split()
255            Ityp = int(s[0])
256            Tminmax = [float(s[1])*1000.,float(s[2])*1000.]
257            Itypes = ['Exponential','Maxwell/Exponential','','Maxwell/Chebyschev','']
258            if Ityp in [1,2,4]:
259                Inst2['Itype'] = Itypes[Ityp-1]
260                Inst2['Tminmax'] = Tminmax
261                Icoeff = []
262                Iesd = []
263                Icovar = []
264                for i in range(3):
265                    s = Iparm['INS  1ICOFF'+str(i+1)].split()
266                    Icoeff += [float(S) for S in s]
267                    s = Iparm['INS  1IECOF'+str(i+1)].split()
268                    Iesd += [float(S) for S in s]
269                NT = 10
270                for i in range(8):
271                    s = Iparm['INS  1IECOR'+str(i+1)]
272                    if i == 7:
273                        NT = 8
274                    Icovar += [float(s[6*j:6*j+6]) for j in range(NT)]
275                Inst2['Icoeff'] = Icoeff
276                Inst2['Iesd'] = Iesd
277                Inst2['Icovar'] = Icovar
278        return [Inst1,Inst2]
279    elif 'E' in DataType:
280        tth = float(Iparm['INS  1 ICONS'])
281        s = Iparm['INS  1PRCF11'].split()
282        names = ['Type','2-theta','XE','YE','ZE','A','B', 'C']
283        codes = [0,0,0,0,0,0,0,0]
284        data.extend([tth,0.0,0.0,0.0,sfloat(s[0]),sfloat(s[1]),sfloat(s[2])])    #A,B,C
285        Iparm1 = makeInstDict(names,data,codes)
286        Iparm1['Bank'] = [Bank,Bank,0]
287        return [Iparm1,{}]
288
289def ReadPowderInstprm(instLines, bank, databanks, rd):
290    '''Read lines from a GSAS-II (new) instrument parameter file
291    similar to G2pwdGUI.OnLoad
292    If instprm file has multiple banks each with header #Bank n: ..., this
293    finds matching bank no. to load - problem with nonmatches?
294   
295    Note that this routine performs a similar role to :meth:`GSASIIdataGUI.GSASII.ReadPowderInstprm`,
296    but that will call a GUI routine for selection when needed. This routine will raise exceptions
297    on errors and will select the first bank when a choice might be appropriate.
298    TODO: refactor to combine the two routines.
299   
300    :param list instLines: strings from GSAS-II parameter file; can be concatenated with ';'
301    :param int  bank: bank number to check when instprm file has '#BANK n:...' strings
302         when bank = n then use parameters; otherwise skip that set. Ignored if BANK n:
303         not present. NB: this kind of instprm file made by a Save all profile command in Instrument Par     ameters
304    :return dict: Inst  instrument parameter dict if OK, or
305             str: Error message if failed
306   
307    (transliterated from GSASIIdataGUI.py:1235 (rev 3008), function of the same name)
308     ''' 
309    if 'GSAS-II' not in instLines[0]:
310        raise ValueError("Not a valid GSAS-II instprm file")
311
312    newItems = []
313    newVals = []
314    Found = False
315    il = 0
316    if bank is None:
317        banklist = set()
318        for S in instLines:
319            if S[0] == '#' and 'Bank' in S:
320                banklist.add(int(S.split(':')[0].split()[1]))
321        # Picks the first bank by default
322        if len(banklist) > 1:
323            bank = sorted(banklist)[0]
324        else:
325            bank = 1
326        rd.powderentry[2] = bank
327    while il < len(instLines):
328        S = instLines[il]
329        if S[0] == '#':
330            if Found:
331                break
332            if 'Bank' in S:
333                if bank == int(S.split(':')[0].split()[1]):
334                    il += 1
335                    S = instLines[il]
336                else:
337                    il += 1
338                    S = instLines[il]
339                    while il < len(instLines) and '#Bank' not in S:
340                        il += 1
341                        if il == len(instLines):
342                            raise ValueError("Bank {} not found in instprm file".format(bank))
343                        S = instLines[il]
344                    continue
345            else:
346                il += 1
347                S = instLines[il]
348        Found = True
349        if '"""' in S:
350            delim = '"""'
351        elif "'''" in S:
352            delim = "'''"
353        else:
354            S = S.replace(' ', '')
355            SS = S.strip().split(';')
356            for s in SS:
357                item, val = s.split(':', 1)
358                newItems.append(item)
359                try:
360                    newVals.append(float(val))
361                except ValueError:
362                    newVals.append(val)
363            il += 1
364            continue
365        # read multiline values, delimited by ''' or """
366        item, val = S.strip().split(':', 1)
367        val = val.replace(delim, '').rstrip()
368        val += '\n'
369        while True:
370            il += 1
371            if il >= len(instLines):
372                break
373            S = instLines[il]
374            if delim in S:
375                val += S.replace(delim, '').rstrip()
376                val += '\n'
377                break
378            else:
379                val += S.rstrip()
380                val += '\n'
381        newItems.append(item)
382        newVals.append(val)
383        il += 1
384    if 'Lam1' in newItems:
385        rd.Sample.update({'Type':'Bragg-Brentano','Shift':[0.,False],'Transparency':[0.,False],
386            'SurfRoughA':[0.,False],'SurfRoughB':[0.,False]})
387    else:
388        rd.Sample.update({'Type':'Debye-Scherrer','Absorption':[0.,False],'DisplaceX':[0.,False],'DisplaceY':[0.,False]})
389    return [makeInstDict(newItems, newVals, len(newVals)*[False]), {}]
390
391def LoadImportRoutines(prefix, errprefix=None, traceback=False):
392    '''Routine to locate GSASII importers matching a prefix string
393    '''
394    if errprefix is None:
395        errprefix = prefix
396
397    readerlist = []
398    import_files = {}
399    if '.' not in sys.path: sys.path.append('.')
400    for path in sys.path:
401        for filename in glob.iglob(os.path.join(path, 'G2'+prefix+'*.py')):
402            pkg = os.path.splitext(os.path.split(filename)[1])[0]
403            if pkg in import_files:
404                G2Print('Warning: file {} hidden by {}'.format(os.path.abspath(filename),import_files[pkg]))
405            else:
406                import_files[pkg] = filename           
407
408    for pkg in sorted(import_files.keys()):
409        try:
410            exec('import '+pkg)
411            for name, value in inspect.getmembers(eval(pkg)):
412                if name.startswith('_'):
413                    continue
414                if inspect.isclass(value):
415                    for method in 'Reader', 'ExtensionValidator', 'ContentsValidator':
416                        if not hasattr(value, method):
417                            break
418                        if not callable(getattr(value, method)):
419                            break
420                    else:
421                        reader = value()
422                        if reader.UseReader:
423                            readerlist.append(reader)
424        except AttributeError:
425            G2Print ('Import_' + errprefix + ': Attribute Error ' + import_files[pkg])
426            if traceback:
427                traceback.print_exc(file=sys.stdout)
428        except Exception as exc:
429            G2Print ('\nImport_' + errprefix + ': Error importing file ' + import_files[pkg])
430            G2Print (u'Error message: {}\n'.format(exc))
431            if traceback:
432                traceback.print_exc(file=sys.stdout)
433    return readerlist
434
435def LoadExportRoutines(parent, traceback=False):
436    '''Routine to locate GSASII exporters
437    '''
438    exporterlist = []
439    import_files = {}
440    if '.' not in sys.path: sys.path.append('.')
441    for path in sys.path:
442        for filename in glob.iglob(os.path.join(path,"G2export*.py")):
443            pkg = os.path.splitext(os.path.split(filename)[1])[0]
444            if pkg in import_files:
445                G2Print('Warning: file {} hidden by {}'
446                                  .format(os.path.abspath(filename),import_files[pkg]))
447            else:
448                import_files[pkg] = filename           
449    # go through the routines and import them, saving objects that
450    # have export routines (method Exporter)
451    for pkg in sorted(import_files.keys()):
452        try:
453            exec('import '+pkg)
454            for clss in inspect.getmembers(eval(pkg)): # find classes defined in package
455                if clss[0].startswith('_'): continue
456                if not inspect.isclass(clss[1]): continue
457                # check if we have the required methods
458                if not hasattr(clss[1],'Exporter'): continue
459                if not callable(getattr(clss[1],'Exporter')): continue
460                if parent is None:
461                    if not hasattr(clss[1],'Writer'): continue
462                else:
463                    if not hasattr(clss[1],'loadParmDict'): continue
464                    if not callable(getattr(clss[1],'loadParmDict')): continue
465                try:
466                    exporter = clss[1](parent) # create an export instance
467                except AttributeError:
468                    pass
469                except Exception as exc:
470                    G2Print ('\nExport init: Error substantiating class ' + clss[0])
471                    G2Print (u'Error message: {}\n'.format(exc))
472                    if traceback:
473                        traceback.print_exc(file=sys.stdout)
474                    continue
475                exporterlist.append(exporter)
476        except AttributeError:
477            G2Print ('Export Attribute Error ' + import_files[pkg])
478            if traceback:
479                traceback.print_exc(file=sys.stdout)
480        except Exception as exc:
481            G2Print ('\nExport init: Error importing file ' + import_files[pkg])
482            G2Print (u'Error message: {}\n'.format(exc))
483            if traceback:
484                traceback.print_exc(file=sys.stdout)
485    return exporterlist
486
487def readColMetadata(imagefile):
488    '''Reads image metadata from a column-oriented metadata table
489    (1-ID style .par file). Called by :func:`GetColumnMetadata`
490   
491    The .par file has any number of columns separated by spaces.
492    The directory for the file must be specified in
493    Config variable :data:`config_example.Column_Metadata_directory`.
494    As an index to the .par file a second "label file" must be specified with the
495    same file root name as the .par file but the extension must be .XXX_lbls (where
496    .XXX is the extension of the image) or if that is not present extension
497    .lbls.
498
499    :param str imagefile: the full name of the image file (with extension, directory optional)
500
501    :returns: a dict with parameter values. Named parameters will have the type based on
502       the specified Python function, named columns will be character strings
503   
504    The contents of the label file will look like this::
505   
506        # define keywords
507        filename:lambda x,y: "{}_{:0>6}".format(x,y)|33,34
508        distance: float | 75
509        wavelength:lambda keV: 12.398425/float(keV)|9
510        pixelSize:lambda x: [74.8, 74.8]|0
511        ISOlikeDate: lambda dow,m,d,t,y:"{}-{}-{}T{} ({})".format(y,m,d,t,dow)|0,1,2,3,4
512        Temperature: float|53
513        FreePrm2: int | 34 | Free Parm2 Label
514        # define other variables
515        0:day
516        1:month
517        2:date
518        3:time
519        4:year
520        7:I_ring
521
522    This file contains three types of lines in any order.
523     * Named parameters are evaluated with user-supplied Python code (see
524       subsequent information). Specific named parameters are used
525       to determine values that are used for image interpretation (see table,
526       below). Any others are copied to the Comments subsection of the Image
527       tree item.
528     * Column labels are defined with a column number (integer) followed by
529       a colon (:) and a label to be assigned to that column. All labeled
530       columns are copied to the Image's Comments subsection.
531     * Comments are any line that does not contain a colon.
532
533    Note that columns are numbered starting at zero.
534
535    Any named parameter may be defined provided it is not a valid integer,
536    but the named parameters in the table have special meanings, as descibed.
537    The parameter name is followed by a colon. After the colon, specify
538    Python code that defines or specifies a function that will be called to
539    generate a value for that parameter.
540
541    Note that several keywords, if defined in the Comments, will be found and
542    placed in the appropriate section of the powder histogram(s)'s Sample
543    Parameters after an integration: ``Temperature``, ``Pressure``, ``Time``,
544    ``FreePrm1``, ``FreePrm2``, ``FreePrm3``, ``Omega``, ``Chi``, and ``Phi``.
545
546    After the Python code, supply a vertical bar (|) and then a list of one
547    more more columns that will be supplied as arguments to that function.
548
549    Note that the labels for the three FreePrm items can be changed by
550    including that label as a third item with an additional vertical bar. Labels
551    will be ignored for any other named parameters.
552   
553    The examples above are discussed here:
554
555    ``filename:lambda x,y: "{}_{:0>6}".format(x,y)|33,34``
556        Here the function to be used is defined with a lambda statement::
557       
558          lambda x,y: "{}_{:0>6}".format(x,y)
559
560        This function will use the format function to create a file name from the
561        contents of columns 33 and 34. The first parameter (x, col. 33) is inserted directly into
562        the file name, followed by a underscore (_), followed by the second parameter (y, col. 34),
563        which will be left-padded with zeros to six characters (format directive ``:0>6``).
564
565        When there will be more than one image generated per line in the .par file, an alternate way to
566        generate list of file names takes into account the number of images generated::
567
568          lambda x,y,z: ["{}_{:0>6}".format(x,int(y)+i) for i in range(int(z))]
569
570        Here a third parameter is used to specify the number of images generated, where
571        the image number is incremented for each image.
572         
573    ``distance: float | 75``
574        Here the contents of column 75 will be converted to a floating point number
575        by calling float on it. Note that the spaces here are ignored.
576       
577    ``wavelength:lambda keV: 12.398425/float(keV)|9``
578        Here we define an algebraic expression to convert an energy in keV to a
579        wavelength and pass the contents of column 9 as that input energy
580       
581    ``pixelSize:lambda x: [74.8, 74.8]|0``
582        In this case the pixel size is a constant (a list of two numbers). The first
583        column is passed as an argument as at least one argument is required, but that
584        value is not used in the expression.
585
586    ``ISOlikeDate: lambda dow,m,d,t,y:"{}-{}-{}T{} ({})".format(y,m,d,t,dow)|0,1,2,3,4``
587        This example defines a parameter that takes items in the first five columns
588        and formats them in a different way. This parameter is not one of the pre-defined
589        parameter names below. Some external code could be used to change the month string
590        (argument ``m``) to a integer from 1 to 12.
591       
592    ``FreePrm2: int | 34 | Free Parm2 Label``
593        In this example, the contents of column 34 will be converted to an integer and
594        placed as the second free-named parameter in the Sample Parameters after an
595        integration. The label for this parameter will be changed to "Free Parm2 Label".
596           
597    **Pre-defined parameter names**
598   
599    =============  =========  ========  =====================================================
600     keyword       required    type      Description
601    =============  =========  ========  =====================================================
602       filename    yes         str or   generates the file name prefix for the matching image
603                               list     file (MyImage001 for file /tmp/MyImage001.tif) or
604                                        a list of file names.
605     polarization  no         float     generates the polarization expected based on the
606                                        monochromator angle, defaults to 0.99.
607       center      no         list of   generates the approximate beam center on the detector
608                              2 floats  in mm, such as [204.8, 204.8].
609       distance    yes        float     generates the distance from the sample to the detector
610                                        in mm
611       pixelSize   no         list of   generates the size of the pixels in microns such as
612                              2 floats  [200.0, 200.0].
613       wavelength  yes        float     generates the wavelength in Angstroms
614    =============  =========  ========  =====================================================
615   
616    '''
617    dir,fil = os.path.split(os.path.abspath(imagefile))
618    imageName,ext = os.path.splitext(fil)
619    if not GSASIIpath.GetConfigValue('Column_Metadata_directory'): return
620    parfiles = glob.glob(os.path.join(GSASIIpath.GetConfigValue('Column_Metadata_directory'),'*.par'))
621    if len(parfiles) == 0:
622        G2Print('Sorry, No Column metadata (.par) file found in '+
623              GSASIIpath.GetConfigValue('Column_Metadata_directory'))
624        return {}
625    for parFil in parfiles: # loop over all .par files (hope just 1) in image dir until image is found
626        parRoot = os.path.splitext(parFil)[0]
627        for e in (ext+'_lbls','.lbls'):
628            if os.path.exists(parRoot+e):
629                lblFil = parRoot+e
630                break
631        else:
632            G2Print('Warning: No labels definitions found for '+parFil)
633            continue
634        labels,lbldict,keyCols,keyExp,errors = readColMetadataLabels(lblFil)
635        if errors:
636            print('Errors in labels file '+lblFil)
637            for i in errors: print('  '+i)
638            continue
639        else:
640            G2Print('Read '+lblFil)
641        # scan through each line in this .par file, looking for the matching image rootname
642        fp = open(parFil,'r')
643        for iline,line in enumerate(fp):
644            items = line.strip().split(' ')
645            nameList = keyExp['filename'](*[items[j] for j in keyCols['filename']])
646            if type(nameList) is str:
647                if nameList != imageName: continue
648                name = nameList
649            else:
650                for name in nameList:
651                    if name == imageName: break # got a match
652                else:
653                    continue
654            # parse the line and finish
655            metadata = evalColMetadataDicts(items,labels,lbldict,keyCols,keyExp)
656            metadata['par file'] = parFil
657            metadata['lbls file'] = lblFil
658            G2Print("Metadata read from {} line {}".format(parFil,iline+1))
659            fp.close()
660            return metadata
661        else:
662            G2Print("Image {} not found in {}".format(imageName,parFil))
663            fp.close()
664            continue
665        fp.close()
666    else:
667        G2Print("Warning: No .par metadata for image {}".format(imageName))
668        return {}
669
670def readColMetadataLabels(lblFil):
671    '''Read the .*lbls file and setup for metadata assignments
672    '''
673    lbldict = {}
674    keyExp = {}
675    keyCols = {}
676    labels = {}
677    errors = []
678    fp = open(lblFil,'r')         # read column labels
679    for iline,line in enumerate(fp): # read label definitions
680        line = line.strip()
681        if not line or line[0] == '#': continue # comments
682        items = line.split(':')
683        if len(items) < 2: continue # lines with no colon are also comments
684        # does this line a definition for a named parameter?
685        key = items[0]
686        try: 
687            int(key)
688        except ValueError: # try as named parameter since not a valid number
689            items = line.split(':',1)[1].split('|')
690            try:
691                f = eval(items[0]) # compile the expression
692                if not callable(f):
693                    errors += ['Expression "{}" for key {} is not a function (line {})'.
694                           format(items[0],key,iline)]
695                    continue
696                keyExp[key] = f
697            except Exception as msg:
698                errors += ['Expression "{}" for key {} is not valid (line {})'.
699                           format(items[0],key,iline)]
700                errors += [str(msg)]
701                continue
702            keyCols[key] = [int(i) for i in items[1].strip().split(',')]
703            if key.lower().startswith('freeprm') and len(items) > 2:
704                labels[key] = items[2]
705            continue
706        if len(items) == 2: # simple column definition
707            lbldict[int(items[0])] = items[1]
708    fp.close()
709    if 'filename' not in keyExp:
710        errors += ["File {} is invalid. No valid filename expression.".format(lblFil)]
711    return labels,lbldict,keyCols,keyExp,errors
712
713def evalColMetadataDicts(items,labels,lbldict,keyCols,keyExp,ShowError=False):
714    '''Evaluate the metadata for a line in the .par file
715    '''
716    metadata = {lbldict[j]:items[j] for j in lbldict}
717    named = {}
718    for key in keyExp:
719        try:
720            res = keyExp[key](*[items[j] for j in keyCols[key]])
721        except:
722            if ShowError:
723                res = "*** error ***"
724            else:
725                continue
726        named[key] = res
727    metadata.update(named)
728    for lbl in labels: # add labels for FreePrm's
729        metadata['label_'+lbl[4:].lower()] = labels[lbl]
730    return metadata
731
732def GetColumnMetadata(reader):
733    '''Add metadata to an image from a column-type metadata file
734    using :func:`readColMetadata`
735   
736    :param reader: a reader object from reading an image
737   
738    '''
739    if not GSASIIpath.GetConfigValue('Column_Metadata_directory'): return
740    parParms = readColMetadata(reader.readfilename)
741    if not parParms: return # check for read failure
742    specialKeys = ('filename',"polarization", "center", "distance", "pixelSize", "wavelength",)
743    reader.Comments = ['Metadata from {} assigned by {}'.format(parParms['par file'],parParms['lbls file'])]
744    for key in parParms:
745        if key in specialKeys+('par file','lbls file'): continue
746        reader.Comments += ["{} = {}".format(key,parParms[key])]
747    if "polarization" in parParms:
748        reader.Data['PolaVal'][0] = parParms["polarization"]
749    else:
750        reader.Data['PolaVal'][0] = 0.99
751    if "center" in parParms:
752        reader.Data['center'] = parParms["center"]
753    if "pixelSize" in parParms:
754        reader.Data['pixelSize'] = parParms["pixelSize"]
755    if "wavelength" in parParms:
756        reader.Data['wavelength'] = parParms['wavelength']
757    else:
758        G2Print('Error: wavelength not defined in {}'.format(parParms['lbls file']))
759    if "distance" in parParms:
760        reader.Data['distance'] = parParms['distance']
761        reader.Data['setdist'] = parParms['distance']
762    else:
763        G2Print('Error: distance not defined in {}'.format(parParms['lbls file']))
764
765def LoadControls(Slines,data):
766    'Read values from a .imctrl (Image Controls) file'
767    cntlList = ['color','wavelength','distance','tilt','invert_x','invert_y','type','Oblique',
768        'fullIntegrate','outChannels','outAzimuths','LRazimuth','IOtth','azmthOff','DetDepth',
769        'calibskip','pixLimit','cutoff','calibdmin','Flat Bkg','varyList','setdist',
770        'PolaVal','SampleAbs','dark image','background image','twoth']
771    save = {}
772    for S in Slines:
773        if S[0] == '#':
774            continue
775        [key,val] = S.strip().split(':',1)
776        if key in ['type','calibrant','binType','SampleShape','color',]:    #strings
777            save[key] = val
778        elif key in ['varyList',]:
779            save[key] = eval(val)   #dictionary
780        elif key in ['rotation']:
781            save[key] = float(val)
782        elif key in ['center',]:
783            if ',' in val:
784                save[key] = eval(val)
785            else:
786                vals = val.strip('[] ').split()
787                save[key] = [float(vals[0]),float(vals[1])] 
788        elif key in cntlList:
789            save[key] = eval(val)
790    data.update(save)
791
792def WriteControls(filename,data):
793    'Write current values to a .imctrl (Image Controls) file'
794    File = open(filename,'w')
795    keys = ['type','color','wavelength','calibrant','distance','center','Oblique',
796            'tilt','rotation','azmthOff','fullIntegrate','LRazimuth','setdist',
797            'IOtth','outChannels','outAzimuths','invert_x','invert_y','DetDepth',
798            'calibskip','pixLimit','cutoff','calibdmin','Flat Bkg','varyList',
799            'binType','SampleShape','PolaVal','SampleAbs','dark image','background image',
800            'twoth']
801    for key in keys:
802        if key not in data:     #uncalibrated!
803            continue
804        File.write(key+':'+str(data[key])+'\n')
805    File.close()
806
807def RereadImageData(ImageReaderlist,imagefile,ImageTag=None,FormatName=''):
808    '''Read a single image with an image importer. This is called to
809    reread an image after it has already been imported, so it is not
810    necessary to reload metadata.
811
812    Based on :func:`GetImageData.GetImageData` which this can replace
813    where imageOnly=True
814
815    :param list ImageReaderlist: list of Reader objects for images
816    :param str imagefile: name of image file
817    :param int/str ImageTag: specifies a particular image to be read from a file.
818      First image is read if None (default).
819    :param str formatName: the image reader formatName
820
821    :returns: an image as a numpy array
822    '''
823    # determine which formats are compatible with this file
824    primaryReaders = []
825    secondaryReaders = []
826    for rd in ImageReaderlist:
827        flag = rd.ExtensionValidator(imagefile)
828        if flag is None: 
829            secondaryReaders.append(rd)
830        elif flag:
831            if not FormatName:
832                primaryReaders.append(rd)
833            elif FormatName == rd.formatName:
834                primaryReaders.append(rd)
835    if len(secondaryReaders) + len(primaryReaders) == 0:
836        G2Print('Error: No matching format for file '+imagefile)
837        raise Exception('No image read')
838    errorReport = ''
839    if not imagefile:
840        return
841    for rd in primaryReaders+secondaryReaders:
842        rd.ReInitialize() # purge anything from a previous read
843        rd.errors = "" # clear out any old errors
844        if not rd.ContentsValidator(imagefile): # rejected on cursory check
845            errorReport += "\n  "+rd.formatName + ' validator error'
846            if rd.errors: 
847                errorReport += ': '+rd.errors
848                continue
849        flag = rd.Reader(imagefile,None,blocknum=ImageTag)
850        if flag: # this read succeeded
851            if rd.Image is None:
852                raise Exception('No image read. Strange!')
853            if GSASIIpath.GetConfigValue('Transpose'):
854                G2Print ('Warning: Transposing Image!')
855                rd.Image = rd.Image.T
856            #rd.readfilename = imagefile
857            return rd.Image
858    else:
859        G2Print('Error reading file '+imagefile)
860        G2Print('Error messages(s)\n'+errorReport)
861        raise Exception('No image read')
862
863def readMasks(filename,masks,ignoreThreshold):
864    '''Read a GSAS-II masks file'''
865    File = open(filename,'r')
866    save = {}
867    oldThreshold = masks['Thresholds'][0]
868    S = File.readline()
869    while S:
870        if S[0] == '#':
871            S = File.readline()
872            continue
873        [key,val] = S.strip().split(':',1)
874        if key in ['Points','Rings','Arcs','Polygons','Frames','Thresholds']:
875            if ignoreThreshold and key == 'Thresholds':
876                S = File.readline() 
877                continue
878            save[key] = eval(val)
879            if key == 'Thresholds':
880                save[key][0] = oldThreshold
881                save[key][1][1] = min(oldThreshold[1],save[key][1][1])
882        S = File.readline()
883    File.close()
884    masks.update(save)
885    # CleanupMasks
886    for key in ['Points','Rings','Arcs','Polygons']:
887        masks[key] = masks.get(key,[])
888        masks[key] = [i for i in masks[key] if len(i)]
889
890def PDFWrite(PDFentry,fileroot,PDFsaves,PDFControls,Inst={},Limits=[]):
891    '''Write PDF-related data (G(r), S(Q),...) into files, as
892    selected.
893
894    :param str PDFentry: name of the PDF entry in the tree. This is
895      used for comments in the file specifying where it came from;
896      it can be arbitrary
897    :param str fileroot: name of file(s) to be written. The extension
898      will be ignored.
899    :param list PDFsaves: flags that determine what type of file will be
900      written:
901      PDFsaves[0], if True writes a I(Q) file with a .iq extension
902      PDFsaves[1], if True writes a S(Q) file with a .sq extension
903      PDFsaves[2], if True writes a F(Q) file with a .fq extension
904      PDFsaves[3], if True writes a G(r) file with a .gr extension
905      PDFsaves[4], if True writes G(r) in a pdfGUI input file with
906      a .gr extension. Note that if PDFsaves[3] and PDFsaves[4] are
907      both True, the pdfGUI overwrites the G(r) file.
908      PDFsaves[5], if True writes F(Q) & g(R) with .fq & .gr extensions
909      overwrites these if selected by option 2, 3 or 4
910    :param dict PDFControls: The PDF parameters and computed results
911    :param dict Inst: Instrument parameters from the PDWR entry used
912      to compute the PDF. Needed only when PDFsaves[4] is True.
913    :param list Limits: Computation limits from the PDWR entry used
914      to compute the PDF. Needed only when PDFsaves[4] is True.
915    '''
916    import scipy.interpolate as scintp
917    fileroot = os.path.splitext(fileroot)[0]
918    if PDFsaves[0]:     #I(Q)
919        iqfilename = fileroot+'.iq'
920        iqdata = PDFControls['I(Q)'][1]
921        iqfxn = scintp.interp1d(iqdata[0],iqdata[1],kind='linear')
922        iqfile = open(iqfilename,'w')
923        iqfile.write('#T I(Q) %s\n'%(PDFentry))
924        iqfile.write('#L Q     I(Q)\n')
925        qnew = np.arange(iqdata[0][0],iqdata[0][-1],0.005)
926        iqnew = zip(qnew,iqfxn(qnew))
927        for q,iq in iqnew:
928            iqfile.write("%15.6g %15.6g\n" % (q,iq))
929        iqfile.close()
930        G2Print (' I(Q) saved to: '+iqfilename)
931
932    if PDFsaves[1]:     #S(Q)
933        sqfilename = fileroot+'.sq'
934        sqdata = PDFControls['S(Q)'][1]
935        sqfxn = scintp.interp1d(sqdata[0],sqdata[1],kind='linear')
936        sqfile = open(sqfilename,'w')
937        sqfile.write('#T S(Q) %s\n'%(PDFentry))
938        sqfile.write('#L Q     S(Q)\n')
939        qnew = np.arange(sqdata[0][0],sqdata[0][-1],0.005)
940        sqnew = zip(qnew,sqfxn(qnew))
941        for q,sq in sqnew:
942            sqfile.write("%15.6g %15.6g\n" % (q,sq))
943        sqfile.close()
944        G2Print (' S(Q) saved to: '+sqfilename)
945
946    if PDFsaves[2]:     #F(Q)
947        fqfilename = fileroot+'.fq'
948        fqdata = PDFControls['F(Q)'][1]
949        fqfxn = scintp.interp1d(fqdata[0],fqdata[1],kind='linear')
950        fqfile = open(fqfilename,'w')
951        fqfile.write('#T F(Q) %s\n'%(PDFentry))
952        fqfile.write('#L Q     F(Q)\n')
953        qnew = np.arange(fqdata[0][0],fqdata[0][-1],0.005)
954        fqnew = zip(qnew,fqfxn(qnew))
955        for q,fq in fqnew:
956            fqfile.write("%15.6g %15.6g\n" % (q,fq))
957        fqfile.close()
958        G2Print (' F(Q) saved to: '+fqfilename)
959
960    if PDFsaves[3]:     #G(R)
961        grfilename = fileroot+'.gr'
962        grdata = PDFControls['G(R)'][1]
963        grfxn = scintp.interp1d(grdata[0],grdata[1],kind='linear')
964        grfile = open(grfilename,'w')
965        grfile.write('#T G(R) %s\n'%(PDFentry))
966        grfile.write('#L R     G(R)\n')
967        rnew = np.arange(grdata[0][0],grdata[0][-1],0.010)
968        grnew = zip(rnew,grfxn(rnew))
969        for r,gr in grnew:
970            grfile.write("%15.6g %15.6g\n" % (r,gr))
971        grfile.close()
972        G2Print (' G(R) saved to: '+grfilename)
973
974    if PDFsaves[4]: #pdfGUI file for G(R)
975        import GSASIImath as G2mth
976        import GSASIIlattice as G2lat       
977        grfilename = fileroot+'.gr'
978        grdata = PDFControls['G(R)'][1]
979        qdata = PDFControls['I(Q)'][1][0]
980        grfxn = scintp.interp1d(grdata[0],grdata[1],kind='linear')
981        grfile = open(grfilename,'w')
982        rnew = np.arange(grdata[0][0],grdata[0][-1],0.010)
983        grnew = zip(rnew,grfxn(rnew))
984
985        grfile.write('[DEFAULT]\n')
986        grfile.write('\n')
987        grfile.write('version = GSAS-II-v'+str(GSASIIpath.GetVersionNumber())+'\n')
988        grfile.write('\n')
989        grfile.write('# input and output specifications\n')
990        grfile.write('dataformat = Qnm\n')
991        grfile.write('inputfile = %s\n'%(PDFControls['Sample']['Name']))
992        grfile.write('backgroundfile = %s\n'%(PDFControls['Sample Bkg.']['Name']))
993        grfile.write('outputtype = gr\n')
994        grfile.write('\n')
995        grfile.write('# PDF calculation setup\n')
996        if 'x' in Inst['Type']:
997            grfile.write('mode = %s\n'%('xray'))
998        elif 'N' in Inst['Type']:
999            grfile.write('mode = %s\n'%('neutron'))
1000        wave = G2mth.getMeanWave(Inst)
1001        grfile.write('wavelength = %.5f\n'%(wave))
1002        formula = ''
1003        for el in PDFControls['ElList']:
1004            formula += el
1005            num = PDFControls['ElList'][el]['FormulaNo']
1006            if num == round(num):
1007                formula += '%d'%(int(num))
1008            else:
1009                formula += '%.2f'%(num)
1010        grfile.write('composition = %s\n'%(formula))
1011        grfile.write('bgscale = %.3f\n'%(-PDFControls['Sample Bkg.']['Mult']))
1012        highQ = 2.*np.pi/G2lat.Pos2dsp(Inst,Limits[1][1])
1013        grfile.write('qmaxinst = %.2f\n'%(highQ))
1014        grfile.write('qmin = %.5f\n'%(qdata[0]))
1015        grfile.write('qmax = %.4f\n'%(qdata[-1]))
1016        grfile.write('rmin = %.2f\n'%(PDFControls['Rmin']))
1017        grfile.write('rmax = %.2f\n'%(PDFControls['Rmax']))
1018        grfile.write('rstep = 0.01\n')
1019        grfile.write('\n')
1020        grfile.write('# End of config '+63*'-')
1021        grfile.write('\n')
1022        grfile.write('#### start data\n')
1023        grfile.write('#S 1\n')
1024        grfile.write('#L r($\\AA$)  G($\\AA^{-2}$)\n')           
1025        for r,gr in grnew:
1026            grfile.write("%15.2F %15.6F\n" % (r,gr))
1027        grfile.close()
1028        G2Print (' G(R) saved to: '+grfilename)
1029       
1030    if len(PDFsaves) > 5 and PDFsaves[5]: #RMCProfile files for F(Q) & g(r) overwrites any above
1031       
1032        fqfilename = fileroot+'.fq'
1033        fqdata = PDFControls['F(Q)'][1]
1034        fqfxn = scintp.interp1d(fqdata[0],fqdata[1],kind='linear')
1035        fqfile = open(fqfilename,'w')
1036        qnew = np.arange(fqdata[0][0],fqdata[0][-1],0.005)
1037        nq = qnew.shape[0]
1038        fqfile.write('%20d\n'%(nq-1))
1039        fqfile.write(fqfilename+'\n')
1040        fqnew = zip(qnew,fqfxn(qnew))
1041        for q,fq in fqnew:
1042            if not q:
1043                continue
1044            fqfile.write("%15.6g %15.6g\n" % (q,fq))
1045        fqfile.close()
1046        G2Print (' F(Q) saved to: '+fqfilename)
1047       
1048        grfilename = fileroot+'.gr'
1049        grdata = PDFControls['g(r)'][1]
1050        grfxn = scintp.interp1d(grdata[0],grdata[1],kind='linear')
1051        grfile = open(grfilename,'w')
1052        rnew = np.arange(grdata[0][0],grdata[0][-1],0.010)
1053        nr = rnew.shape[0]
1054        grfile.write('%20d\n'%(nr-1))
1055        grfile.write(grfilename+'\n')
1056        grnew = zip(rnew,grfxn(rnew))
1057        for r,gr in grnew:
1058            if not r:
1059                continue
1060            grfile.write("%15.6g %15.6g\n" % (r,gr))
1061        grfile.close()
1062        G2Print (' G(R) saved to: '+grfilename)
1063
Note: See TracBrowser for help on using the repository browser.