source: trunk/GSASIIIO.py @ 239

Last change on this file since 239 was 239, checked in by vondreele, 11 years ago

add 'any image file' to image file menu
add calibration skip & dmin to image data dictionary
fix to ellipse fitting
fix Pilatus reading - OK for 100K, not sure for 2M
now a image sizexy - 2 items for x & y sizes

  • Property svn:keywords set to Date Author Revision URL Id
File size: 44.4 KB
Line 
1"""GSASIIIO: functions for IO of data
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2011-01-13 19:34:07 +0000 (Thu, 13 Jan 2011) $
6# $Author: vondreele $
7# $Revision: 239 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 239 2011-01-13 19:34:07Z vondreele $
10########### SVN repository information ###################
11import wx
12import math
13import numpy as np
14import cPickle
15import sys
16import random as ran
17import GSASIIpath
18import GSASIIgrid as G2gd
19import GSASIIspc as G2spc
20import GSASIIlattice as G2lat
21import GSASIIElem as G2el
22import os.path as ospath
23
24def sfloat(S):
25    if S.strip():
26        return float(S)
27    else:
28        return 0.0
29
30def sint(S):
31    if S.strip():
32        return int(S)
33    else:
34        return 0
35
36def SelectPowderData(self, filename):
37    """Selects banks of data from a filename of any GSAS powder data format
38    Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
39    Returns - a list of banks to be read; each entry in list is a tuple containing:
40    filename: same as input filename
41    Pos: position for start of data; record just after BANK record
42    Bank: the BANK record
43    """
44    File = open(filename,'Ur')
45    Title = '''
46First line of this file:
47'''+File.readline()
48    dlg = wx.MessageDialog(self, Title, 'Is this the file you want?', 
49    wx.YES_NO | wx.ICON_QUESTION)
50    try:
51        result = dlg.ShowModal()
52    finally:
53        dlg.Destroy()
54    if result == wx.ID_NO: return (0,0)
55    Temperature = 300
56       
57    self.IparmName = GetInstrumentFile(self,filename)
58    if self.IparmName:
59        Iparm = GetInstrumentData(self.IparmName)
60    else:
61        Iparm = {}                                               #Assume CuKa lab data if no iparm file
62        Iparm['INS   HTYPE '] = 'PXC '
63        Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
64        Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
65        Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
66        Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
67    S = 1
68    Banks = []
69    Pos = []
70    FoundData = []
71    Comments = []
72    wx.BeginBusyCursor()
73    try:
74        while S:
75            S = File.readline()
76            if S[:1] != '#':
77                if S[:4] == 'BANK':
78                    Banks.append(S)
79                    Pos.append(File.tell())
80            else:
81                Comments.append(S[:-1])
82                if 'Temp' in S:
83                    Temperature = float(S[:-1].split()[-1])
84        File.close()
85    finally:
86        wx.EndBusyCursor()
87    if Comments:
88       print 'Comments on file:'
89       for Comment in Comments: print Comment
90    if Banks:
91        result = [0]
92        if len(Banks) >= 2:
93            dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
94            try:
95                if dlg.ShowModal() == wx.ID_OK:
96                    result = dlg.GetSelections()
97                else:
98                    result = []
99            finally:
100                dlg.Destroy()
101        for i in result:
102            FoundData.append((filename,Pos[i],Banks[i]))
103    else:
104        dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
105        try:
106            result = dlg.ShowModal()
107        finally:
108            dlg.Destroy()
109    return FoundData,Iparm,Comments,Temperature
110
111def GetInstrumentFile(self,filename):
112    import os.path as op
113    dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', wx.OPEN)
114    if self.dirname: 
115        dlg.SetDirectory(self.dirname)
116        Tname = filename[:filename.index('.')]+'.prm'
117        if op.exists(Tname):
118            self.IparmName = Tname       
119    if self.IparmName: dlg.SetFilename(self.IparmName)
120    filename = ''
121    try:
122        if dlg.ShowModal() == wx.ID_OK:
123            filename = dlg.GetPath()
124    finally:
125        dlg.Destroy()
126    return filename
127
128def GetInstrumentData(IparmName):
129    file = open(IparmName, 'Ur')
130    S = 1
131    Iparm = {}
132    while S:
133        S = file.readline()
134        Iparm[S[:12]] = S[12:-1]
135    return Iparm
136   
137def GetPowderPeaks(fileName):
138    sind = lambda x: math.sin(x*math.pi/180.)
139    asind = lambda x: 180.*math.asin(x)/math.pi
140    Cuka = 1.54052
141    File = open(fileName,'Ur')
142    Comments = []
143    peaks = []
144    S = File.readline()
145    while S:
146        if S[:1] == '#':
147            Comments.append(S[:-1])
148        else:
149            item = S.split()
150            if len(item) == 1:
151                peaks.append([float(item[0]),1.0])
152            elif len(item) > 1:
153                peaks.append([float(item[0]),float(item[0])])
154        S = File.readline()
155    File.close()
156    if Comments:
157       print 'Comments on file:'
158       for Comment in Comments: print Comment
159    Peaks = []
160    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
161        for peak in peaks:
162            dsp = peak[0]
163            sth = Cuka/(2.0*dsp)
164            if sth < 1.0:
165                tth = 2.0*asind(sth)
166            else:
167                break
168            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
169    else:                                   #2-thetas - assume Cuka (for now)
170        for peak in peaks:
171            tth = peak[0]
172            dsp = Cuka/(2.0*sind(tth/2.0))
173            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
174    return Comments,Peaks
175
176def GetPawleyPeaks(filename):
177    rt2ln2x2 = 2.35482
178    File = open(filename,'Ur')
179    PawleyPeaks = []
180    S = File.readline()         #skip header
181    S = File.readline()
182    item = S.split()
183    while S:
184        h,k,l = int(item[0]),int(item[1]),int(item[2])
185        mult = int(item[3])
186        tth = float(item[5])
187        sig = float(item[6])/rt2ln2x2
188        Iobs = float(item[7])*mult
189        PawleyPeaks.append([h,k,l,mult,tth,sig,False,Iobs,0.0,[]])
190        S = File.readline()
191        item = S.split()
192        if item[3] == '-100.0000':       #find trailer
193            break
194    File.close()
195    return PawleyPeaks
196   
197def GetHKLData(filename):
198    print 'Reading: '+filename
199    File = open(filename,'Ur')
200    HKLref = []
201    HKLmin = [1000,1000,1000]
202    HKLmax = [0,0,0]
203    FoMax = 0
204    ifFc = False
205    S = File.readline()
206    while '#' in S[0]:        #get past comments if any
207        S = File.readline()       
208    if '_' in S:         #cif style .hkl file
209        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
210            S = File.readline()       
211        S = File.readline()             #get past 'loop_' line
212        pos = 0
213        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
214        while S:
215            if '_' in S:
216                if 'index_h' in S:
217                    hpos = pos
218                elif 'index_k' in S:
219                    kpos = pos
220                elif 'index_l' in S:
221                    lpos = pos
222                elif 'F_squared_meas' in S:
223                    Fosqpos = pos
224                elif 'F_squared_calc' in S:
225                    Fcsqpos = pos
226                elif 'F_squared_sigma' in S:
227                    sigpos = pos
228                pos += 1
229            else:
230                data = S.split()
231                if data:                    #avoid blank lines
232                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
233                    h,k,l = HKL
234                    Fosq = float(data[Fosqpos])
235                    if sigpos != -1:
236                        sigFosq = float(data[sigpos])
237                    else:
238                        sigFosq = 1.
239                    if Fcsqpos != -1:
240                        Fcsq = float(data[Fcsqpos])
241                        if Fcsq:
242                            ifFc = True
243                    else:
244                        Fcsq = 0.
245                       
246                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
247                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
248                    FoMax = max(FoMax,Fosq)
249                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fcp, Fcpp & phase
250            S = File.readline()
251    else:                   #dumb h,k,l,Fo,sigFo .hkl file
252        while S:
253            h,k,l,Fo,sigFo = S.split()
254            HKL = np.array([int(h),int(k),int(l)])
255            h,k,l = HKL
256            Fo = float(Fo)
257            sigFo = float(sigFo)
258            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
259            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
260            FoMax = max(FoMax,Fo)
261            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
262            S = File.readline()
263    File.close()
264    return HKLref,HKLmin,HKLmax,FoMax,ifFc
265
266def GetPowderData(filename,Pos,Bank,DataType):
267    '''Reads one BANK of data from GSAS raw powder data file
268    input:
269    filename: GSAS raw powder file dataname
270    Pos: start of data in file just after BANK record
271    Bank: the BANK record
272    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
273    returns: list [x,y,e,yc,yb]
274    x: np.array of x-axis values
275    y: np.array of powder pattern intensities
276    w: np.array of w=sig(intensity)^2 values
277    yc: np.array of calc. intensities (zero)
278    yb: np.array of calc. background (zero)
279    yd: np.array of obs-calc profiles
280    '''
281    print 'Reading: '+filename
282    print 'Bank:    '+Bank[:-1]
283    if 'FXYE' in Bank:
284        return GetFXYEdata(filename,Pos,Bank,DataType)
285    elif 'FXY' in Bank:
286        return GetFXYdata(filename,Pos,Bank,DataType)
287    elif 'ESD' in Bank:
288        return GetESDdata(filename,Pos,Bank,DataType)
289    elif 'STD' in Bank:
290        return GetSTDdata(filename,Pos,Bank,DataType)
291    else:
292        return GetSTDdata(filename,Pos,Bank,DataType)
293    return []
294
295def GetFXYEdata(filename,Pos,Bank,DataType):
296    File = open(filename,'Ur')
297    File.seek(Pos)
298    x = []
299    y = []
300    w = []
301    S = File.readline()
302    while S and S[:4] != 'BANK':
303        vals = S.split()
304        if DataType[2] == 'C':
305            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
306        elif DataType[2] == 'T':
307            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
308        f = float(vals[1])
309        if f <= 0.0:
310            y.append(0.0)
311            w.append(1.0)
312        else:
313            y.append(float(vals[1]))
314            w.append(1.0/float(vals[2])**2)
315        S = File.readline()
316    File.close()
317    N = len(x)
318    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
319   
320def GetFXYdata(filename,Pos,Bank,DataType):
321    File = open(filename,'Ur')
322    File.seek(Pos)
323    x = []
324    y = []
325    w = []
326    S = File.readline()
327    while S and S[:4] != 'BANK':
328        vals = S.split()
329        if DataType[2] == 'C':
330            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
331        elif DataType[2] == 'T':
332            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
333        f = float(vals[1])
334        if f > 0.0:
335            y.append(f)
336            w.append(1.0/f)
337        else:             
338            y.append(0.0)
339            w.append(1.0)
340        S = File.readline()
341    File.close()
342    N = len(x)
343    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
344   
345def GetESDdata(filename,Pos,Bank,DataType):
346    File = open(filename,'Ur')
347    cons = Bank.split()
348    if DataType[2] == 'C':
349        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
350        step = float(cons[6])/100.0
351    elif DataType[2] == 'T':
352        start = float(cons[5])/1000.0              #TOF: from musec to millisec
353        step = float(cons[6])/1000.0
354    File.seek(Pos)
355    x = []
356    y = []
357    w = []
358    S = File.readline()
359    j = 0
360    while S and S[:4] != 'BANK':
361        for i in range(0,80,16):
362            xi = start+step*j
363            yi = sfloat(S[i:i+8])
364            ei = sfloat(S[i+8:i+16])
365            x.append(xi)
366            if yi > 0.0:
367                y.append(yi)
368                w.append(1.0/ei**2)
369            else:             
370                y.append(0.0)
371                w.append(1.0)
372            j += 1
373        S = File.readline()
374    File.close()
375    N = len(x)
376    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
377
378def GetSTDdata(filename,Pos,Bank,DataType):
379    File = open(filename,'Ur')
380    cons = Bank.split()
381    Nch = cons[2]
382    if DataType[2] == 'C':
383        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
384        step = float(cons[6])/100.0
385    elif DataType[2] == 'T':
386        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
387        step = float(cons[6])/1000.0
388    File.seek(Pos)
389    x = []
390    y = []
391    w = []
392    S = File.readline()
393    j = 0
394    while S and S[:4] != 'BANK':
395        for i in range(0,80,8):
396            xi = start+step*j
397            ni = max(sint(S[i:i+2]),1)
398            yi = max(sfloat(S[i+2:i+8]),0.0)
399            if yi:
400                ei = math.sqrt(yi*ni)
401            else:
402                yi = 0.0
403                ei = 1.0
404            j += 1
405            if j < Nch:
406                x.append(xi)
407                y.append(yi)
408                w.append(1.0/ei**2)
409        S = File.readline()
410    File.close()
411    N = len(x)
412    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
413   
414def CheckImageFile(self,imagefile):
415    if not ospath.exists(imagefile):
416        dlg = wx.FileDialog(self, 'Bad image file name; choose name', '.', '',\
417        'MAR345 (*.mar3450;*.mar2300)|*.mar3450;*.mar2300|ADSC Image (*.img)\
418        |*.img|Detector tif (*.tif;*.tiff)|*.tif;*.tiff|GE Image sum (*.sum)\
419        |*.sum|GE Image avg (*.avg)|*.avg|GE Image raw (*)|*|All files (*.*)|*.*',wx.OPEN)
420        if self.dirname:
421            dlg.SetDirectory(self.dirname)
422        try:
423            dlg.SetFilename(ospath.split(imagefile)[1])
424            if dlg.ShowModal() == wx.ID_OK:
425                self.dirname = dlg.GetDirectory()
426                imagefile = dlg.GetPath()
427            else:
428                imagefile = False
429        finally:
430            dlg.Destroy()
431    return imagefile
432       
433def GetImageData(self,imagefile,imageOnly=False):       
434    ext = ospath.splitext(imagefile)[1]
435    Comments = []
436    if ext == '.tif':
437        Comments,Data,Size,Image = GetTifData(imagefile)
438    elif ext == '.img':
439        Comments,Data,Size,Image = GetImgData(imagefile)
440        Image[0][0] = 0
441    elif ext == '.mar3450' or ext == '.mar2300':
442        Comments,Data,Size,Image = GetMAR345Data(imagefile)
443    elif ext in ['.sum','.avg','']:
444        Comments,Data,Size,Image = GetGEsumData(imagefile)
445    elif ext == '.G2img':
446        return GetG2Image(imagefile)
447    if imageOnly:
448        return Image
449    else:
450        return Comments,Data,Size,Image
451       
452def PutG2Image(filename,image):
453    File = open(filename,'wb')
454    cPickle.dump(image,File,1)
455    File.close()
456    return
457   
458def GetG2Image(filename):
459    File = open(filename,'rb')
460    image = cPickle.load(File)
461    File.close()
462    return image
463   
464def GetGEsumData(filename,imageOnly=False):
465    import array as ar
466    if not imageOnly:
467        print 'Read GE sum file: ',filename   
468    File = open(filename,'rb')
469    size = 2048
470    row = 0
471    pos = 0
472    if '.sum' in filename:
473        head = ['GE detector sum data from APS 1-ID',]
474    if '.avg' in filename:
475        head = ['GE detector avg data from APS 1-ID',]
476    else:
477        head = ['GE detector raw data from APS 1-ID',]
478        pos = 8192
479    image = np.zeros(shape=(size,size),dtype=np.int32)
480    while row < size:
481        File.seek(pos)
482        if '.sum' in filename:
483            line = ar.array('f',File.read(4*size))
484            pos += 4*size
485        elif '.avg' in filename:
486            line = ar.array('H',File.read(2*size))
487            pos += 2*size
488        else:
489            line = ar.array('H',File.read(2*size))
490            pos += 2*size
491        image[row] = np.asarray(line)
492        row += 1
493    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':[size,size]} 
494    File.close()   
495    if imageOnly:
496        return image
497    else:
498        return head,data,size,image
499       
500def GetImgData(filename,imageOnly=False):
501    import struct as st
502    import array as ar
503    if not imageOnly:
504        print 'Read ADSC img file: ',filename
505    File = open(filename,'rb')
506    head = File.read(511)
507    lines = head.split('\n')
508    head = []
509    center = [0,0]
510    for line in lines[1:-2]:
511        line = line.strip()[:-1]
512        if line:
513            if 'SIZE1' in line:
514                size = int(line.split('=')[1])
515            elif 'WAVELENGTH' in line:
516                wave = float(line.split('=')[1])
517            elif 'BIN' in line:
518                if line.split('=')[1] == '2x2':
519                    pixel=(102,102)
520                else:
521                    pixel = (51,51)
522            elif 'DISTANCE' in line:
523                distance = float(line.split('=')[1])
524            elif 'CENTER_X' in line:
525                center[0] = float(line.split('=')[1])
526            elif 'CENTER_Y' in line:
527                center[1] = float(line.split('=')[1])
528            head.append(line)
529    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
530    image = []
531    row = 0
532    pos = 512
533    image = np.zeros(shape=(size,size),dtype=np.int32)   
534    while row < size:
535        File.seek(pos)
536        line = ar.array('H',File.read(2*size))
537        image[row] = np.asarray(line)
538        row += 1
539        pos += 2*size
540    File.close()
541    if imageOnly:
542        return image
543    else:
544        return lines[1:-2],data,size,image
545       
546def GetMAR345Data(filename,imageOnly=False):
547    import array as ar
548    import struct as st
549    try:
550        import pack_f as pf
551    except:
552        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
553                               caption="Import Error",
554                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
555        msg.ShowModal()
556        return None,None,None,None
557
558    if not imageOnly:
559        print 'Read Mar345 file: ',filename
560    File = open(filename,'rb')
561    head = File.read(4095)
562    numbers = st.unpack('<iiiiiiiiii',head[:40])
563    lines = head[128:].split('\n')
564    head = []
565    for line in lines:
566        line = line.strip()
567        if 'PIXEL' in line:
568            values = line.split()
569            pixel = (int(values[2]),int(values[4]))     #in microns
570        elif 'WAVELENGTH' in line:
571            wave = float(line.split()[1])
572        elif 'DISTANCE' in line:
573            distance = float(line.split()[1])           #in mm
574        elif 'CENTER' in line:
575            values = line.split()
576            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
577        if line: 
578            head.append(line)
579    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
580    for line in head:
581        if 'FORMAT' in line[0:6]:
582            items = line.split()
583            size = int(items[1])
584    pos = 4096
585    data['size'] = [size,size]
586    File.seek(pos)
587    line = File.read(8)
588    while 'CCP4' not in line:       #get past overflow list for now
589        line = File.read(8)
590        pos += 8
591    pos += 37
592    File.seek(pos)
593    raw = File.read()
594    File.close()
595    image = np.zeros(shape=(size,size),dtype=np.int32)
596    image = pf.pack_f(len(raw),raw,size,image)
597    if imageOnly:
598        return image.T              #transpose to get it right way around
599    else:
600        return head,data,size,image.T
601   
602def GetTifData(filename,imageOnly=False):
603    # only works for APS Perkin-Elmer or MAR detector data files in "TIFF" format that are readable by Fit2D
604    import struct as st
605    import array as ar
606    File = open(filename,'rb')
607    dataType = 5
608    try:
609        Meta = open(filename+'.metadata','Ur')
610        head = Meta.readlines()
611        for line in head:
612            line = line.strip()
613            if 'dataType' in line:
614                dataType = int(line.split('=')[1])
615        Meta.close()
616    except IOError:
617        print 'no metadata file found - will try to read file anyway'
618        head = 'no metadata file found'
619    tag = File.read(3)
620    if tag != 'II*':
621        lines = ['not a detector tiff file',]
622        return lines,0,0
623    origSize,Ityp = st.unpack('<ii',File.read(8))
624    File.seek(30)
625    finalSize = st.unpack('<i',File.read(4))[0]
626    if finalSize:
627        sizeScale = origSize/finalSize
628    else:
629        sizeScale = 1
630        finalSize = origSize
631    if Ityp == 0:
632        File.seek(62)
633        S = File.read(32)
634        if 'PILATUS' in S:
635            dataType = 0
636            tifType = 'Pilatus'
637            if '2M' in S:
638                sizexy = [1475,1679]
639            elif '100K' in S:
640                sizexy = [487,195]
641            pixy = (172,172)
642            pos = 4096
643            if not imageOnly:
644                print 'Read Pilatus tiff file: ',filename
645        else:
646            sizexy = [1536,1536]
647            sizeScale = 1
648            tifType = 'Gold'
649            pixy = (150,150)
650            pos = 64
651            if not imageOnly:
652                print 'Read Gold tiff file:',filename
653    elif Ityp == 1:
654        tifType = 'PE'
655        sizexy = [finalSize,finalSize]
656        pixy = (200*sizeScale,200*sizeScale)
657        pos = 8
658        if not imageOnly:
659            print 'Read APS PE-detector tiff file: ',filename
660    elif Ityp == 3328:
661        tifType = 'MAR'
662        sizexy = [finalSize,finalSize]
663        pixy = (79*sizeScale,79*sizeScale)
664        pos = 4096
665        if not imageOnly:
666            print 'Read MAR CCD tiff file: ',filename
667    else:
668        lines = 'unknown tif type'
669        return lines,0,0
670    image = np.zeros(shape=(sizexy[1],sizexy[0]),dtype=np.int32)
671    row = 0
672    while row < sizexy[1]:
673        File.seek(pos)
674        if 'PE' in tifType or 'Pilatus' in tifType: 
675            if dataType == 5:
676                line = ar.array('f',File.read(4*sizexy[0]))
677            else:
678                line = ar.array('L',File.read(4*sizexy[0]))
679            pos += 4*sizexy[0]
680        elif 'MAR' in tifType or 'Gold' in tifType:
681            line = ar.array('H',File.read(2*sizexy[0]))
682            pos += 2*sizexy[0]
683        image[row] = np.asarray(line)
684        row += 1
685    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
686    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
687    File.close()   
688    if imageOnly:
689        return image
690    else:
691        return head,data,sizexy[0],image
692   
693def ProjFileOpen(self):
694    file = open(self.GSASprojectfile,'rb')
695    print 'load from file: ',self.GSASprojectfile
696    wx.BeginBusyCursor()
697    try:
698        while True:
699            try:
700                data = cPickle.load(file)
701            except EOFError:
702                break
703            datum = data[0]
704            print 'load: ',datum[0]
705           
706            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
707            if 'PWDR' in datum[0]:               
708                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
709            else:
710                self.PatternTree.SetItemPyData(Id,datum[1])
711            for datus in data[1:]:
712                print '    load: ',datus[0]
713                               
714                sub = self.PatternTree.AppendItem(Id,datus[0])
715                self.PatternTree.SetItemPyData(sub,datus[1])
716                                               
717            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
718                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
719                if Data['setDefault']:
720                    self.imageDefault = Data               
721        file.close()
722       
723    finally:
724        wx.EndBusyCursor()
725    print 'project load successful'
726    self.NewPlot = True
727   
728def ProjFileSave(self):
729    if not self.PatternTree.IsEmpty():
730        file = open(self.GSASprojectfile,'wb')
731        print 'save to file: ',self.GSASprojectfile
732        wx.BeginBusyCursor()
733        try:
734            item, cookie = self.PatternTree.GetFirstChild(self.root)
735            while item:
736                data = []
737                name = self.PatternTree.GetItemText(item)
738                print 'save: ',name
739                data.append([name,self.PatternTree.GetItemPyData(item)])
740                item2, cookie2 = self.PatternTree.GetFirstChild(item)
741                while item2:
742                    name = self.PatternTree.GetItemText(item2)
743                    print '    save: ',name
744                    data.append([name,self.PatternTree.GetItemPyData(item2)])
745                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
746                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
747                cPickle.dump(data,file,1)
748            file.close()
749        finally:
750            wx.EndBusyCursor()
751        print 'project save successful'
752       
753def SaveIntegration(self,PickId,data):
754    azms = self.Integrate[1]
755    X = self.Integrate[2][:-1]
756    Xminmax = [X[0],X[-1]]
757    N = len(X)
758    Id = self.PatternTree.GetItemParent(PickId)
759    name = self.PatternTree.GetItemText(Id)
760    name = name.replace('IMG ','PWDR ')
761    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
762    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
763    codes = [0 for i in range(11)]
764    parms = ['PXC',data['wavelength'],0.0,0.0,1.0,-1.0,0.3,0.0,1.0,0.0,0.0]
765    Azms = [(azms[i+1]+azms[i])/2. for i in range(len(azms)-1)]
766    for i,azm in enumerate(Azms):
767        item, cookie = self.PatternTree.GetFirstChild(self.root)
768        Id = 0
769        while item:
770            Name = self.PatternTree.GetItemText(item)
771            if name == Name:
772                Id = item
773            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
774        parms[10] = azm
775        Y = self.Integrate[0][i]
776        W = np.sqrt(Y)
777        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],'DisplaceX':[0.0,False],
778            'DisplaceY':[0.0,False],'Diffuse':[],'Temperature':300.,'Pressure':1.0,'Humidity':0.0,'Voltage':0.0,'Force':0.0}
779        if Id:
780            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
781            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
782            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
783            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms,codes,names])
784            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
785            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
786            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
787        else:
788            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(azm))
789            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
790            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
791            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
792            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms,codes,names])
793            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
794            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
795            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
796            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
797        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
798    self.PatternTree.SelectItem(Id)
799    self.PatternTree.Expand(Id)
800    self.PatternId = Id
801           
802def powderFxyeSave(self,powderfile):
803    file = open(powderfile,'w')
804    prm = open(powderfile.strip('fxye')+'prm','w')      #old style GSAS parm file
805    print 'save powder pattern to file: ',powderfile
806    wx.BeginBusyCursor()
807    Inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
808                    self.PickId, 'Instrument Parameters'))[1]
809    if len(Inst) == 11:             #single wavelength
810        lam1 = Inst[1]
811        lam2 = 0.0
812        GU,GV,GW = Inst[4:7]
813        LX,LY = Inst[7:9]
814        SL = HL = Inst[9]/2.0   
815    else:                           #Ka1 & Ka2
816        lam1 = Inst[1]
817        lam2 = Inst[2]
818        GU,GV,GW = Inst[6:9]
819        LX,LY = Inst[9:11]
820        SL = HL = Inst[11]/2.0   
821    prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
822    prm.write( 'INS   BANK      1                                                               '+'\n')
823    prm.write( 'INS   HTYPE   PXCR                                                              '+'\n')
824    prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(lam1,lam2))
825    prm.write( 'INS  1 IRAD     0                                                               '+'\n')
826    prm.write( 'INS  1I HEAD                                                                    '+'\n')
827    prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
828    prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
829    prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(GU,GV,GW,0.0))
830    prm.write(('INS  1PRCF12     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(LX,LY,SL,HL))
831    prm.close()
832    try:
833        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
834        file.write(powderfile+'\n')
835        file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
836            100.*x[0],100.*(x[1]-x[0])))
837        s = list(np.sqrt(1./np.array(w)))       
838        XYW = zip(x,y,s)
839        for X,Y,S in XYW:
840            file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,S))
841        file.close()
842    finally:
843        wx.EndBusyCursor()
844    print 'powder pattern file written'
845       
846def powderXyeSave(self,powderfile):
847    file = open(powderfile,'w')
848    print 'save powder pattern to file: ',powderfile
849    wx.BeginBusyCursor()
850    try:
851        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
852        XYW = zip(x,y,w)
853        for X,Y,W in XYW:
854            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
855        file.close()
856    finally:
857        wx.EndBusyCursor()
858    print 'powder pattern file written'
859   
860def PeakListSave(self,file,peaks):
861    print 'save peak list to file: ',self.peaklistfile
862    if not peaks:
863        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
864        try:
865            result = dlg.ShowModal()
866        finally:
867            dlg.Destroy()
868        return
869    for peak in peaks:
870        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
871            (peak[0],peak[2],peak[4],peak[6]))
872    print 'peak list saved'
873             
874def IndexPeakListSave(self,peaks):
875    file = open(self.peaklistfile,'wa')
876    print 'save index peak list to file: ',self.peaklistfile
877    wx.BeginBusyCursor()
878    try:
879        if not peaks:
880            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
881            try:
882                result = dlg.ShowModal()
883            finally:
884                dlg.Destroy()
885            return
886        for peak in peaks:
887            file.write("%12.6f\n" % (peak[7]))
888        file.close()
889    finally:
890        wx.EndBusyCursor()
891    print 'index peak list saved'
892   
893def ReadEXPPhase(self,filename):
894    file = open(filename, 'Ur')
895    Phase = {}
896    S = 1
897    Expr = [{},{},{},{},{},{},{},{},{}]
898    while S:
899        S = file.readline()
900        if 'EXPR NPHAS' in S[:12]:
901            Num = S[12:-1].count('0')
902            NPhas = S[12:-1].split()
903        if 'CRS' in S[:3]:
904            N = int(S[3:4])-1
905            Expr[N][S[:12]] = S[12:-1]
906    file.close()
907    PNames = []
908    for n,N in enumerate(NPhas):
909        if N != '0':
910            result = n
911            key = 'CRS'+str(n+1)+'    PNAM'
912            PNames.append(Expr[n][key])
913    if Num < 8:
914        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
915        try:
916            if dlg.ShowModal() == wx.ID_OK:
917                result = dlg.GetSelection()
918        finally:
919            dlg.Destroy()       
920    EXPphase = Expr[result]
921    keyList = EXPphase.keys()
922    keyList.sort()
923    SGData = {}
924    if NPhas[result] == '1':
925        Ptype = 'nuclear'
926    elif NPhas[result] in ['2','3']:
927        Ptype = 'magnetic'
928    elif NPhas[result] == '4':
929        Ptype = 'macromolecular'
930    elif NPhas[result] == '10':
931        Ptype = 'Pawley'
932    for key in keyList:
933        if 'PNAM' in key:
934           PhaseName = EXPphase[key].strip()
935        elif 'ABC   ' in key:
936            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
937        elif 'ANGLES' in key:
938            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
939        elif 'SG SYM' in key:
940            SpGrp = EXPphase[key][:15].strip()
941            E,SGData = G2spc.SpcGroup(SpGrp)
942    Atoms = []
943    if Ptype == 'nuclear':
944        for key in keyList:
945            if 'AT' in key:
946                if key[11:] == 'A':
947                    S = EXPphase[key]
948                elif key[11:] == 'B':
949                    S += EXPphase[key]
950                    Atom = [S[50:58].strip(),S[:10].strip(),'',
951                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
952                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
953                    if Atom[9] == 'I':
954                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
955                    elif Atom[9] == 'A':
956                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
957                            float(S[88:98]),float(S[98:108]),
958                            float(S[108:118]),float(S[118:128])]
959                    XYZ = Atom[3:6]
960                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
961                    Atom.append(ran.randint(0,sys.maxint))
962                    Atoms.append(Atom)
963    elif Ptype == 'macromolecular':
964        for key in keyList:
965            if 'AT' in key[6:8]:
966                S = EXPphase[key]
967                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
968                    S[46:51].strip(),S[:8].strip(),'',
969                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
970                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
971                XYZ = Atom[6:9]
972                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
973                Atom.append(ran.randint(0,sys.maxint))
974                Atoms.append(Atom)
975    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
976    Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
977        'Cell':[False,]+abc+angles+[Volume,]}
978    Phase['Atoms'] = Atoms
979    Phase['Drawing'] = {}
980    Phase['Histograms'] = {}
981    return Phase
982       
983def ReadPDBPhase(filename):
984    EightPiSq = 8.*math.pi**2
985    file = open(filename, 'Ur')
986    Phase = {}
987    Title = ''
988    Compnd = ''
989    Atoms = []
990    A = np.zeros(shape=(3,3))
991    S = file.readline()
992    while S:
993        Atom = []
994        if 'TITLE' in S[:5]:
995            Title = S[10:72].strip()
996            S = file.readline()
997        elif 'COMPND    ' in S[:10]:
998            Compnd = S[10:72].strip()
999            S = file.readline()
1000        elif 'CRYST' in S[:5]:
1001            abc = S[7:34].split()
1002            angles = S[34:55].split()
1003            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1004                float(angles[0]),float(angles[1]),float(angles[2])]
1005            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1006            AA,AB = G2lat.cell2AB(cell)
1007            SpGrp = S[55:65]
1008            E,SGData = G2spc.SpcGroup(SpGrp)
1009            if E: 
1010                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1011                print ' N.B.: make sure spaces separate axial fields in symbol' 
1012                print G2spc.SGErrors(E)
1013                return None
1014            SGlines = G2spc.SGPrint(SGData)
1015            for line in SGlines: print line
1016            S = file.readline()
1017        elif 'SCALE' in S[:5]:
1018            V = (S[10:41].split())
1019            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1020            S = file.readline()
1021        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1022            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1023            XYZ = np.inner(AB,XYZ)
1024            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1025            Uiso = float(S[61:67])/EightPiSq
1026            Type = S[12:14].upper()
1027            if Type[0] in '123456789':
1028                Type = Type[1:]
1029            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1030                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1031                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1032            S = file.readline()
1033            if 'ANISOU' in S[:6]:
1034                Uij = S[30:72].split()
1035                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1036                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1037                Atom = Atom[:14]+Uij
1038                Atom[12] = 'A'
1039                S = file.readline()
1040            Atom.append(ran.randint(0,sys.maxint))
1041            Atoms.append(Atom)
1042        else:           
1043            S = file.readline()
1044    file.close()
1045    if Title:
1046        PhaseName = Title
1047    elif Compnd:
1048        PhaseName = Compnd
1049    else:
1050        PhaseName = 'None'
1051    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1052        'Cell':[False,]+cell+[Volume,]}
1053    Phase['Atoms'] = Atoms
1054    Phase['Drawing'] = {}
1055    Phase['Histograms'] = {}
1056   
1057    return Phase
1058   
1059def ReadCIFPhase(filename):
1060    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1061    file = open(filename, 'Ur')
1062    Phase = {}
1063    Title = ospath.split(filename)[-1]
1064    print '\n Reading cif file: ',Title
1065    Compnd = ''
1066    Atoms = []
1067    A = np.zeros(shape=(3,3))
1068    S = file.readline()
1069    while S:
1070        if '_symmetry_space_group_name_H-M' in S:
1071            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1072            E,SGData = G2spc.SpcGroup(SpGrp)
1073            if E:
1074                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1075                print ' N.B.: make sure spaces separate axial fields in symbol' 
1076                print G2spc.SGErrors(E)
1077                return None
1078            S = file.readline()
1079        elif '_cell' in S:
1080            if '_cell_length_a' in S:
1081                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1082            elif '_cell_length_b' in S:
1083                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1084            elif '_cell_length_c' in S:
1085                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1086            elif '_cell_angle_alpha' in S:
1087                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1088            elif '_cell_angle_beta' in S:
1089                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1090            elif '_cell_angle_gamma' in S:
1091                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1092            S = file.readline()
1093        elif 'loop_' in S:
1094            labels = {}
1095            i = 0
1096            while S:
1097                S = file.readline()
1098                if '_atom_site' in S.strip()[:10]:
1099                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1100                    i += 1
1101                else:
1102                    break
1103            if labels:
1104                if 'aniso_label' not in labels:
1105                    while S:
1106                        atom = ['','','',0,0,0,0,'','','I',0.01,0,0,0,0,0,0]
1107                        S.strip()
1108                        if len(S.split()) != len(labels):
1109                            if 'loop_' in S:
1110                                break
1111                            S += file.readline().strip()
1112                        data = S.split()
1113                        if len(data) != len(labels):
1114                            break
1115                        for key in labels:
1116                            if key == 'type_symbol':
1117                                atom[1] = data[labels[key]]
1118                            elif key == 'label':
1119                                atom[0] = data[labels[key]]
1120                            elif key == 'fract_x':
1121                                atom[3] = float(data[labels[key]].split('(')[0])
1122                            elif key == 'fract_y':
1123                                atom[4] = float(data[labels[key]].split('(')[0])
1124                            elif key == 'fract_z':
1125                                atom[5] = float(data[labels[key]].split('(')[0])
1126                            elif key == 'occupancy':
1127                                atom[6] = float(data[labels[key]].split('(')[0])
1128                            elif key == 'thermal_displace_type':
1129                                if data[labels[key]].lower() == 'uiso':
1130                                    atom[9] = 'I'
1131                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1132                                else:
1133                                    atom[9] = 'A'
1134                                    atom[10] = 0.0
1135                                   
1136                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1137                        atom.append(ran.randint(0,sys.maxint))
1138                        Atoms.append(atom)
1139                        S = file.readline()
1140                else:
1141                    while S:
1142                        S.strip()
1143                        data = S.split()
1144                        if len(data) != len(labels):
1145                            break
1146                        name = data[labels['aniso_label']]
1147                        for atom in Atoms:
1148                            if name == atom[0]:
1149                                for i,uname in enumerate(anisoNames):
1150                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1151                        S = file.readline()
1152                                                                       
1153        else:           
1154            S = file.readline()
1155    file.close()
1156    if Title:
1157        PhaseName = Title
1158    else:
1159        PhaseName = 'None'
1160    SGlines = G2spc.SGPrint(SGData)
1161    for line in SGlines: print line
1162    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1163    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1164    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1165        'Cell':[False,]+cell+[Volume,]}
1166    Phase['Atoms'] = Atoms
1167    Phase['Drawing'] = {}
1168    Phase['Histograms'] = {}
1169   
1170    return Phase
Note: See TracBrowser for help on using the repository browser.