source: trunk/imports/G2pwd_CIF.py @ 1069

Last change on this file since 1069 was 1069, checked in by toby, 8 years ago

finally a working CIF exporter (export mechanism needs configuring); work around windows file names in PyCifRW & missing unicode char

  • Property svn:eol-style set to native
File size: 13.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2012-02-13 11:33:35 -0600 (Mon, 13 Feb 2012) $
4# $Author: vondreele & toby $
5# $Revision: 482 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/G2importphase.py $
7# $Id: G2importphase.py 482 2012-02-13 17:33:35Z vondreele $
8########### SVN repository information ###################
9# routines to read in structure factors from a CIF
10#
11import sys
12import numpy as np
13import os.path
14import GSASIIIO as G2IO
15import CifFile as cif # PyCifRW from James Hester
16import urllib
17import GSASIIpath
18GSASIIpath.SetVersionNumber("$Revision: 810 $")
19
20class CIFpwdReader(G2IO.ImportPowderData):
21    'Routines to import powder data from a CIF file'
22    def __init__(self):
23        super(self.__class__,self).__init__( # fancy way to self-reference
24            extensionlist=('.CIF','.cif'),
25            strictExtension=False,
26            formatName = 'CIF',
27            longFormatName = 'Powder data from CIF'
28            )
29    # Validate the contents
30    def ContentsValidator(self, filepointer):
31        for i,line in enumerate(filepointer):
32            if i >= 1000: break
33            ''' Encountered only blank lines or comments in first 1000
34            lines. This is unlikely, but assume it is CIF since we are
35            even less likely to find a file with nothing but hashes and
36            blank lines'''
37            line = line.strip()
38            if len(line) == 0:
39                continue # ignore blank lines
40            elif line.startswith('#'):
41                continue # ignore comments
42            elif line.startswith('data_'):
43                return True
44            else:
45                return False # found something else
46        return True
47
48    def Reader(self,filename,filepointer, ParentFrame=None, **kwarg):
49        # Define lists of data names used for holding powder diffraction data
50        xDataItems = (  # "x-axis" data names"
51            ('_pd_meas_2theta_range_min', '_pd_meas_2theta_range_max', '_pd_meas_2theta_range_inc'),
52            ('_pd_proc_2theta_range_min', '_pd_proc_2theta_range_max', '_pd_proc_2theta_range_inc'),
53            '_pd_meas_2theta_scan',
54            '_pd_meas_time_of_flight',
55            '_pd_proc_2theta_corrected',
56            #'_pd_proc_d_spacing',
57            #'_pd_proc_energy_incident',
58            #'_pd_proc_energy_detection',
59            #'_pd_proc_recip_len_q',
60            #'_pd_proc_wavelength',
61        )
62        intDataItems = ( # intensity axis data names
63            '_pd_meas_counts_total',
64            #'_pd_meas_counts_background',
65            #'_pd_meas_counts_container',
66            '_pd_meas_intensity_total',
67            #'_pd_meas_intensity_background',
68            #'_pd_meas_intensity_container',
69            '_pd_proc_intensity_net',
70            '_pd_proc_intensity_total',
71            #'_pd_proc_intensity_bkg_calc',
72            #'_pd_proc_intensity_bkg_fix',
73            '_pd_calc_intensity_net',   # allow computed patterns as input data?
74            '_pd_calc_intensity_total',
75            )
76
77        ESDDataItems = ( # items that can be used to compute uncertainties for obs data
78            '_pd_proc_ls_weight',
79            '_pd_meas_counts_total'
80            )
81       
82        ModDataItems = ( # items that modify the use of the data
83            '_pd_meas_step_count_time',
84            '_pd_meas_counts_monitor',
85            '_pd_meas_intensity_monitor',
86            '_pd_proc_intensity_norm',
87            '_pd_proc_intensity_incident',
88            )
89        rdbuffer = kwarg.get('buffer')
90        cf = None
91        choicelist = None
92        selections = None
93        # reload previously saved values
94        if self.repeat and rdbuffer is not None:
95            cf = rdbuffer.get('lastcif')
96            choicelist = rdbuffer.get('choicelist')
97            print 'debug: Reuse previously parsed CIF'
98            selections = rdbuffer.get('selections')
99        try:
100            if cf is None:
101                self.ShowBusy() # this can take a while
102                ciffile = 'file:'+urllib.pathname2url(filename)
103                #cf = cif.ReadCif(ciffile)
104                fp = open(ciffile,'r')             # patch: open file to avoid windows bug
105                cf = cif.ReadCif(fp)
106                fp.close()
107                self.DoneBusy()
108        except Exception as detail:
109            print self.formatName+' read error:'+str(detail) # for testing
110            import traceback
111            traceback.print_exc(file=sys.stdout)
112            return False
113               
114        # scan all blocks for sets of data
115        if choicelist is None:
116            choicelist = []
117            for blk in cf.keys():
118                blkkeys = [k.lower() for k in cf[blk].keys()] # save a list of the data items, since we will use it often
119                # scan through block for x items
120                xldict = {}
121                for x in xDataItems:
122                    if type(x) is tuple: # check for the presence of all three items that define a range of data
123                        if not all([i in blkkeys for i in x]): continue
124                        try:
125                            items = [float(cf[blk][xi]) for xi in x]
126                            l = 1 + int(0.5 + (items[1]-items[0])/items[2])
127                        except:
128                            continue
129                    else:
130                        if x not in blkkeys: continue
131                        l = len(cf[blk][x])
132                    if xldict.get(l) is None:
133                        xldict[l] = [x]
134                    else:
135                        xldict[l].append(x)
136                # now look for matching intensity items
137                yldict = {}
138                suldict = {}
139                for y in intDataItems:
140                    if y in blkkeys:
141                        l = len(cf[blk][y])
142                        if yldict.get(l) is None:
143                            yldict[l] = [y]
144                        else:
145                            yldict[l].append(y)
146                        # now check if the first item has an uncertainty
147                        if cif.get_number_with_esd(cf[blk][y][0])[1] is None: continue
148                        if suldict.get(l) is None:
149                            suldict[l] = [y]
150                        else:
151                            suldict[l].append(y)
152                for y in ESDDataItems:
153                    if y in blkkeys:
154                        l = len(cf[blk][y])
155                        if suldict.get(l) is None:
156                            suldict[l] = [y]
157                        else:
158                            suldict[l].append(y)
159                modldict = {}
160                for y in ModDataItems:
161                    if y in blkkeys:
162                        l = len(cf[blk][y])
163                        if modldict.get(l) is None:
164                            modldict[l] = [y]
165                        else:
166                            modldict[l].append(y)
167                for l in xldict:
168                    if yldict.get(l) is None: continue
169                    choicelist.append([blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])])
170                    #print blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])
171        if not choicelist:
172            selblk = None # no block to choose
173            return False
174        elif len(choicelist) == 1: # only one choice
175            selblk = 0
176        elif self.repeat and selections is not None:
177            # we were called to repeat the read
178            print 'debug: repeat #',self.repeatcount,'selection',selections[self.repeatcount]
179            selblk = selections[self.repeatcount]
180            self.repeatcount += 1
181            if self.repeatcount >= len(selections): self.repeat = False
182        else:                       # choose from options
183            # compile a list of choices for the user
184            choices = []
185            for blk,l,x,y,su,mod in choicelist:
186                sx = x[0]
187                if len(x) > 1: sx += '...'
188                sy = y[0]
189                if len(y) > 1: sy += '...'
190                choices.append(
191                    'Block '+str(blk)+', '+str(l)+' points. X='+sx+' & Y='+sy
192                    )
193            selections = self.MultipleBlockSelector(
194                choices,
195                ParentFrame=ParentFrame,
196                title='Select dataset(s) to read from the list below',
197                size=(600,100),
198                header='Dataset Selector')
199            if len(selections) == 0: return False
200            selblk = selections[0] # select first in list
201            if len(selections) > 1: # prepare to loop through again
202                self.repeat = True
203                self.repeatcount = 1
204                if rdbuffer is not None:
205                    rdbuffer['selections'] = selections
206                    rdbuffer['lastcif'] = cf # save the parsed cif for the next loop
207                    rdbuffer['choicelist'] = choicelist # save the parsed choices for the future
208
209        # got a selection, now read it
210        blk,l,xch,ych,such,modch = choicelist[selblk]
211        xi,yi,sui,modi = 0,0,0,0
212        if len(xch) > 1 or len(ych) > 1 or len(such) > 1 or len(modch) > 0:
213            choices = []
214            chlbls = []
215            chlbls.append('Select the scanned data item')
216            xchlst = []
217            for x in xch:
218                if type(x) is tuple:
219                    xchlst.append(x[0])
220                else:
221                    xchlst.append(x)
222            choices.append(xchlst)
223            chlbls.append('Select the intensity data item')
224            choices.append(ych)
225            chlbls.append('Select the data item to be used for weighting')
226            choices.append(such)
227            chlbls.append('Divide intensities by data item')
228            choices.append(['none']+modch)
229            res = self.MultipleChoicesDialog(choices,chlbls)
230            if not res: return False
231            xi,yi,sui,modi = res
232
233        # now read in the values
234        # x-values
235        xcf = xch[xi]
236        if type(xcf) is tuple:
237            vals = [float(cf[blk][xi]) for xi in xcf]
238            x = np.array(
239                [(i*vals[2] + vals[0]) for i in range(1 + int(0.5 + (vals[1]-vals[0])/vals[2]))]
240                )
241        else:
242            vl = []
243            for val in cf[blk].get(xcf,'?'):
244                v,e = cif.get_number_with_esd(val)
245                if v is None: # not parsed
246                    vl.append(np.NaN)
247                else:
248                    vl.append(v)
249            x = np.array(vl)
250        # y-values
251        ycf = ych[yi]
252        vl = []
253        for val in cf[blk].get(ycf,'?'):
254            v,e = cif.get_number_with_esd(val)
255            if v is None: # not parsed
256                vl.append(np.NaN)
257            else:
258                vl.append(v)
259        y = np.array(vl)
260        # weights
261        if sui == -1:
262            # no weights
263            vl = np.zeros(len(x)) + 1.
264        else:
265            sucf = such[sui]
266            if sucf ==  '_pd_proc_ls_weight':
267                for val in cf[blk].get(sucf,'?'):
268                    v,e = cif.get_number_with_esd(val)
269                    if v is None: # not parsed
270                        vl.append(0.)
271                    else:
272                        vl.append(v)
273            elif sucf ==  '_pd_meas_counts_total':
274                for val in cf[blk].get(sucf,'?'):
275                    v,e = cif.get_number_with_esd(val)
276                    if v is None: # not parsed
277                        vl.append(0.)
278                    elif v <= 0:
279                        vl.append(1.)
280                    else:
281                        vl.append(1./v)
282            else:
283                for val in cf[blk].get(sucf,'?'):
284                    v,e = cif.get_number_with_esd(val)
285                    if v is None: # not parsed
286                        vl.append(0.0)
287                    elif v <= 0:
288                        vl.append(0.0)
289                    else:
290                        vl.append(1./(v*v))
291        w = np.array(vl)
292        # intensity modification factor
293        if modi >= 1:
294            ycf = modch[modi-1]
295            vl = []
296            for val in cf[blk].get(ycf,'?'):
297                v,e = cif.get_number_with_esd(val)
298                if v is None: # not parsed
299                    vl.append(np.NaN)
300                else:
301                    vl.append(v)
302            y /= np.array(vl)
303            w /= np.array(vl)
304        N = len(x)
305        self.powderdata = [
306                np.array(x), # x-axis values
307                np.array(y), # powder pattern intensities
308                np.array(w), # 1/sig(intensity)^2 values (weights)
309                np.zeros(N), # calc. intensities (zero)
310                np.zeros(N), # calc. background (zero)
311                np.zeros(N), # obs-calc profiles
312            ]
313        self.powderentry[0] = filename
314        #self.powderentry[1] = pos # bank offset (N/A here)
315        self.powderentry[2] = 1 # xye file only has one bank
316        self.idstring = os.path.basename(filename) + ': ' + blk
317        if cf[blk].get('_diffrn_radiation_probe'):
318            if cf[blk]['_diffrn_radiation_probe'] == 'neutron':
319                self.instdict['type'] = 'SNC'
320                #if cf[blk].get('_pd_meas_time_of_flight'): self.instdict['type'] = 'SNT' # not supported yet
321            else:
322                self.instdict['type'] = 'SXC'
323        if cf[blk].get('_diffrn_radiation_wavelength'):
324            val = cf[blk]['_diffrn_radiation_wavelength']
325            wl = []
326            if type(val) is list:
327                for v in val:
328                    w,e = cif.get_number_with_esd(v)
329                    if w: wl.append(w)
330            else:
331                w,e = cif.get_number_with_esd(val)
332                if w: wl.append(w)
333            if wl: self.instdict['wave'] = wl
334        if cf[blk].get('_diffrn_ambient_temperature'):
335            val = cf[blk]['_diffrn_ambient_temperature']
336            w,e = cif.get_number_with_esd(val)
337            if w: 
338                self.Sample['Temperature'] = w
339        return True
Note: See TracBrowser for help on using the repository browser.