source: trunk/GSASIIfiles.py @ 4018

Last change on this file since 4018 was 4018, checked in by vondreele, 2 years ago

set default to Bragg-Brentano for dual wavelength source & using GSAS-II insparm file on import of powder data - now works same as for old GSAS iparm files

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