source: trunk/imports/G2pwd_CIF.py @ 3786

Last change on this file since 3786 was 3786, checked in by vondreele, 3 years ago

fix problem of import ElementTable? inside spyder
allow import o q-steped powder data from a cif file
fix problem of indexing after load incommensurate phase in Unit cells

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 14.2 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2019-01-17 20:31:32 +0000 (Thu, 17 Jan 2019) $
4# $Author: vondreele $
5# $Revision: 3786 $
6# $URL: trunk/imports/G2pwd_CIF.py $
7# $Id: G2pwd_CIF.py 3786 2019-01-17 20:31:32Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2pwd_CIF: CIF powder data*
11------------------------------------
12
13Routine to read in powder data from a CIF.
14
15'''
16from __future__ import division, print_function
17import numpy as np
18import os.path
19import GSASIIobj as G2obj
20import GSASIIIO as G2IO
21import CifFile as cif # PyCifRW from James Hester
22import GSASIIpath
23asind = lambda x: 180.*np.arcsin(x)/np.pi
24GSASIIpath.SetVersionNumber("$Revision: 3786 $")
25
26class CIFpwdReader(G2obj.ImportPowderData):
27    'Routines to import powder data from a CIF file'
28    def __init__(self):
29        super(self.__class__,self).__init__( # fancy way to self-reference
30            extensionlist=('.CIF','.cif'),
31            strictExtension=False,
32            formatName = 'CIF',
33            longFormatName = 'Powder data from CIF'
34            )
35    # Validate the contents
36    def ContentsValidator(self, filename):
37        'Use standard CIF validator'
38        fp = open(filename,'r')
39        return self.CIFValidator(fp)
40        fp.close()
41
42    def Reader(self,filename, ParentFrame=None, **kwarg):
43        '''Read powder data from a CIF.
44        If multiple datasets are requested, use self.repeat and buffer caching.
45        '''
46
47        # Define lists of data names used for holding powder diffraction data
48        # entries of a type that are not implemented are commented out in case
49        # we will want them later.
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        if cf is None:
100            if GSASIIpath.GetConfigValue('debug'): print("Starting parse of {} as CIF".format(filename))
101            cf = G2obj.ReadCIF(filename)
102            if GSASIIpath.GetConfigValue('debug'): print ("CIF file parsed")
103        # scan all blocks for sets of data
104        if choicelist is None:
105            choicelist = []
106            for blk in cf.keys():
107                blkkeys = [k.lower() for k in cf[blk].keys()] # save a list of the data items, since we will use it often
108                # scan through block for x items
109                xldict = {}
110                for x in xDataItems:
111                    if type(x) is tuple: # check for the presence of all three items that define a range of data
112                        if not all([i in blkkeys for i in x]): continue
113                        try:
114                            items = [float(cf[blk][xi]) for xi in x]
115                            l = 1 + int(0.5 + (items[1]-items[0])/items[2])
116                        except:
117                            continue
118                    else:
119                        if x not in blkkeys: continue
120                        l = len(cf[blk][x])
121                    if xldict.get(l) is None:
122                        xldict[l] = [x]
123                    else:
124                        xldict[l].append(x)
125                # now look for matching intensity items
126                yldict = {}
127                suldict = {}
128                for y in intDataItems:
129                    if y in blkkeys:
130                        l = len(cf[blk][y])
131                        if yldict.get(l) is None:
132                            yldict[l] = [y]
133                        else:
134                            yldict[l].append(y)
135                        # now check if the first item has an uncertainty
136                        if cif.get_number_with_esd(cf[blk][y][0])[1] is None: continue
137                        if suldict.get(l) is None:
138                            suldict[l] = [y]
139                        else:
140                            suldict[l].append(y)
141                for y in ESDDataItems:
142                    if y in blkkeys:
143                        l = len(cf[blk][y])
144                        if suldict.get(l) is None:
145                            suldict[l] = [y]
146                        else:
147                            suldict[l].append(y)
148                modldict = {}
149                for y in ModDataItems:
150                    if y in blkkeys:
151                        l = len(cf[blk][y])
152                        if modldict.get(l) is None:
153                            modldict[l] = [y]
154                        else:
155                            modldict[l].append(y)
156                for l in xldict:
157                    if yldict.get(l) is None: continue
158                    choicelist.append([blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])])
159                    #print blk,l,xldict[l],yldict[l],suldict.get(l,[]),modldict.get(l,[])
160            print ("CIF file scanned for blocks with data")
161        if not choicelist:
162            selblk = None # no block to choose
163            self.errors = "No powder diffraction blocks found"
164            return False
165        elif len(choicelist) == 1: # only one choice
166            selblk = 0
167        elif self.repeat and selections is not None:
168            # we were called to repeat the read
169            print ('debug: repeat #',self.repeatcount,'selection',selections[self.repeatcount])
170            selblk = selections[self.repeatcount]
171            self.repeatcount += 1
172            if self.repeatcount >= len(selections): self.repeat = False
173        else:                       # choose from options
174            # compile a list of choices for the user
175            choices = []
176            for blk,l,x,y,su,mod in choicelist:
177                sx = x[0]
178                if len(x) > 1: sx += '...'
179                sy = y[0]
180                if len(y) > 1: sy += '...'
181                choices.append(
182                    'Block '+str(blk)+', '+str(l)+' points. X='+sx+' & Y='+sy
183                    )
184            selections = G2IO.MultipleBlockSelector(
185                choices,
186                ParentFrame=ParentFrame,
187                title='Select dataset(s) to read from the list below',
188                size=(600,100),
189                header='Dataset Selector')
190            if len(selections) == 0:
191                self.errors = "Abort: block not selected"
192                return False
193            selblk = selections[0] # select first in list
194            if len(selections) > 1: # prepare to loop through again
195                self.repeat = True
196                self.repeatcount = 1
197                if rdbuffer is not None:
198                    rdbuffer['selections'] = selections
199                    rdbuffer['lastcif'] = cf # save the parsed cif for the next loop
200                    rdbuffer['choicelist'] = choicelist # save the parsed choices for the future
201
202        # got a selection, now read it
203        # do we need to ask which fields to read?
204        blk,l,xch,ych,such,modch = choicelist[selblk]
205        xi,yi,sui,modi = 0,0,0,0
206        if len(xch) > 1 or len(ych) > 1 or len(such) > 1 or len(modch) > 0:
207            choices = []
208            chlbls = []
209            chlbls.append('Select the scanned data item')
210            xchlst = []
211            for x in xch:
212                if type(x) is tuple:
213                    xchlst.append(x[0])
214                else:
215                    xchlst.append(x)
216            choices.append(xchlst)
217            chlbls.append('Select the intensity data item')
218            choices.append(ych)
219            chlbls.append('Select the data item to be used for weighting')
220            choices.append(such)
221            chlbls.append('Divide intensities by data item')
222            choices.append(['none']+modch)
223            res = G2IO.MultipleChoicesSelector(choices,chlbls)
224            if not res:
225                self.errors = "Abort: data items not selected"
226                return False
227            xi,yi,sui,modi = res
228
229            # now read in the values
230            # x-values
231            self.powderentry[0] = filename
232            #self.powderentry[1] = pos # bank offset (N/A here)
233            #self.powderentry[2] = 1 # xye file only has one bank
234            self.idstring = os.path.basename(filename) + ': ' + blk
235            if cf[blk].get('_diffrn_radiation_probe'):
236                if cf[blk]['_diffrn_radiation_probe'] == 'neutron':
237                    self.instdict['type'] = 'PNC'
238                    #if cf[blk].get('_pd_meas_time_of_flight'): self.instdict['type'] = 'PNT' # not supported yet
239                else:
240                    self.instdict['type'] = 'PXC'
241            if cf[blk].get('_diffrn_radiation_wavelength'):
242                val = cf[blk]['_diffrn_radiation_wavelength']
243                wl = []
244                if type(val) is list:
245                    for v in val:
246                        w,e = cif.get_number_with_esd(v)
247                        if w: wl.append(w)
248                else:
249                    w,e = cif.get_number_with_esd(val)
250                    if w: wl.append(w)
251                if wl: 
252                    if len(wl) > 1:
253                        self.instdict['wave'] = wl
254                    else:
255                        self.instdict['wave'] = wl[0]
256            if cf[blk].get('_diffrn_ambient_temperature'):
257                val = cf[blk]['_diffrn_ambient_temperature']
258                w,e = cif.get_number_with_esd(val)
259                if w: 
260                    self.Sample['Temperature'] = w
261        xcf = xch[xi]
262        if type(xcf) is tuple:
263            vals = [float(cf[blk][ixi]) for ixi in xcf]
264            x = np.array([(i*vals[2] + vals[0]) for i in range(1 + int(0.5 + (vals[1]-vals[0])/vals[2]))])
265        else:
266            vl = []
267            for val in cf[blk].get(xcf,'?'):
268                v,e = cif.get_number_with_esd(val)
269                if v is None: # not parsed
270                    vl.append(np.NaN)
271                else:
272                    vl.append(v)
273            x = np.array(vl)
274            if 'recip_len_q' in xcf and 'wave' in self.instdict:
275                wl = self.instdict['wave']
276                x = 2.*asind(wl*x/(4.*np.pi))
277        # y-values
278        ycf = ych[yi]
279        vl = []
280        v2 = []
281        for val in cf[blk].get(ycf,'?'):
282            v,e = cif.get_number_with_esd(val)
283            if v is None: # not parsed
284                vl.append(np.NaN)
285                v2.append(np.NaN)
286            else:
287                vl.append(v)
288                if e is None:
289                    v2.append(np.sqrt(v))
290                else:
291                    v2.append(max(e,1.0))
292        y = np.array(vl)
293        w = 1./np.array(v2)**2
294        # weights
295        if sui == -1:
296            # no weights
297            vl = np.zeros(len(x)) + 1.
298        else:
299            vl = []
300            sucf = such[sui]
301            if sucf ==  '_pd_proc_ls_weight':
302                for val in cf[blk].get(sucf,'?'):
303                    v,e = cif.get_number_with_esd(val)
304                    if v is None: # not parsed
305                        vl.append(0.)
306                    else:
307                        vl.append(v)
308            elif sucf ==  '_pd_meas_counts_total':
309                for val in cf[blk].get(sucf,'?'):
310                    v,e = cif.get_number_with_esd(val)
311                    if v is None: # not parsed
312                        vl.append(0.)
313                    elif v <= 0:
314                        vl.append(1.)
315                    else:
316                        vl.append(1./v)
317            else:
318                for val in cf[blk].get(sucf,'?'):
319                    v,e = cif.get_number_with_esd(val)
320                    if v is None or e is None: # not parsed or no ESD
321                        vl.append(np.NaN)
322                    elif e <= 0:
323                        vl.append(1.)
324                    else:
325                        vl.append(1./(e*e))
326#            w = np.array(vl)
327        # intensity modification factor
328        if modi >= 1:
329            ycf = modch[modi-1]
330            vl = []
331            for val in cf[blk].get(ycf,'?'):
332                v,e = cif.get_number_with_esd(val)
333                if v is None: # not parsed
334                    vl.append(np.NaN)
335                else:
336                    vl.append(v)
337            y /= np.array(vl)
338            w /= np.array(vl)
339        N = len(x)
340        print ("CIF file, read from selected block")
341
342        self.errors = "Error while storing read values"
343        self.powderdata = [
344                np.array(x), # x-axis values
345                np.array(y), # powder pattern intensities
346                np.array(w), # 1/sig(intensity)^2 values (weights)
347                np.zeros(N), # calc. intensities (zero)
348                np.zeros(N), # calc. background (zero)
349                np.zeros(N), # obs-calc profiles
350            ]
351        return True
Note: See TracBrowser for help on using the repository browser.