source: TI02-CSML/trunk/csml/csmllibs/csmlcrs.py @ 2665

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI02-CSML/trunk/csml/csmllibs/csmlcrs.py@2665
Revision 2665, 11.9 KB checked in by domlowe, 13 years ago (diff)

used dictionary to store fileextract info that was getting read from disk each time

Line 
1import string
2import sys         
3
4class CRSystem(object):
5    def __init__(self, srsName, axes):
6        self.srsName=srsName
7        self.axes=axes
8        self.srsDimension=len(self.axes)
9        self.axisLabels=''
10        for axis in self.axes:
11            self.axisLabels=self.axisLabels + axis + ' '
12        self.axisLabels=self.axisLabels[:-1]
13        self.timeAxis=None
14        self.lonAxis=None
15        self.latAxis=None
16        self.units=None
17   
18    def getAxisLabels(self):
19        labelList=[]
20        for item in self.axisLabels.split():
21            labelList.append(item)
22        return labelList
23
24class CRSCatalogue(object):
25    def __init__(self):
26        # dictionary to hold known CRSystems:
27        self.systems={}
28        self.pseudosystems={}
29       
30        # define WGS84 CRS:       
31        crs=CRSystem(srsName='WGS84', axes =['time', 'latitude','longitude'])
32        crs.lonAxis=2
33        crs.latAxis=1
34        crs.timeAxis=0
35        self.systems['WGS84']=crs
36       
37       
38         #time only crs
39        crs=CRSystem(srsName='ndg:crs:t', axes =['time'])
40        crs.timeAxis=0
41        crs.lonAxis=None
42        crs.latAxis=None
43        self.systems['ndg:crs:t']=crs
44       
45        # define lon lat pressure time CRS:       
46        crs=CRSystem(srsName='ndg:crs:xypt', axes =['longitude', 'latitude','pressure','time'])
47        crs.lonAxis=0
48        crs.latAxis=1
49        crs.timeAxis=3
50        self.systems['ndg:crs:xypt']=crs
51       
52        # define lon lat height time CRS:       
53        crs=CRSystem(srsName='ndg:crs:xyht', axes =['longitude', 'latitude','height','time'])
54        crs.lonAxis=0
55        crs.latAxis=1
56        crs.timeAxis=3
57        self.systems['ndg:crs:xyht']=crs
58       
59        # define lon lat time CRS:
60        crs=CRSystem(srsName='ndg:crs:xyt', axes =['longitude', 'latitude','time'])
61        self.systems['ndg:crs:xyt']=crs
62        crs.lonAxis=0
63        crs.latAxis=1
64        crs.timeAxis=2
65       
66        # define lon lat height CRS:
67        crs=CRSystem(srsName='ndg:crs:xyh', axes =['longitude', 'latitude','height'])
68        self.systems['ndg:crs:xyh']=crs
69        crs.lonAxis=0
70        crs.latAxis=1
71        crs.timeAxis=None
72       
73       
74        # define height time CRS:
75        crs=CRSystem(srsName='ndg:crs:ht', axes =['height', 'time'])
76        self.systems['ndg:crs:ht']=crs
77        crs.lonAxis=None
78        crs.latAxis=None
79        crs.timeAxis=1
80       
81        #define unknown 2D CRS:
82        crs=CRSystem(srsName='ndg:crs:undefined2d', axes=['tba','time'])
83        self.pseudosystems['ndg:crs:undefined2d']=crs
84        crs.timeAxis=1
85        crs.lonAxis=None
86        crs.latAxis=None
87       
88        #define unknown 3D CRS:
89        crs=CRSystem(srsName='ndg:crs:undefinedll3d', axes=['longitude','latitude','tba'])
90        self.pseudosystems['ndg:crs:undefinedll3d']=crs
91        crs.lonAxis=0
92        crs.latAxis=1
93        crs.timeAxis=None
94       
95        #define unknown 4D CRS:
96        crs=CRSystem(srsName='ndg:crs:undefinedllt4d', axes=['longitude','latitude','tba','time'])
97        self.pseudosystems['ndg:crs:undefinedllt4d']=crs
98        crs.lonAxis=0
99        crs.latAxis=1
100        crs.timeAxis=3
101       
102        #define unknown 1D CRS:
103        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
104        self.systems['ndg:crs:unknown']=crs
105                       
106        #define unknown 3D CRS:
107        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
108        self.systems['ndg:crs:unknown2d']=crs
109       
110        #define unknown 3D CRS:
111        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
112        self.systems['ndg:crs:unknown3d']=crs
113       
114        #define unknown 4D CRS:
115        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
116        self.systems['ndg:crs:unknown4d']=crs
117       
118           
119    def getCRS(self, srsName, labels=None):
120        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
121        #labels are only required if it is a ndg:pseudoCRS
122        try:
123            return self.systems[srsName]
124        except:
125            if srsName=='ndg:pseudoCRS':
126                labellist=[]
127                try:
128                    for label in labels.split():
129                        labellist.append(label)
130                except:
131                        labellist=labels
132            #need to find the right pseudosystem.
133            #at the moment a simple test on number of dimensions works, but need to make this more generic.
134            for psystem in self.pseudosystems:
135                try:
136                    if len(self.pseudosystems[psystem].axes) == len (labellist):
137                        return self.pseudosystems[psystem]
138                except:
139                    pass
140            else:
141                return self.systems['ndg:crs:unknown']
142   
143    def getUnitType(self, unit):
144        #unittype='unknown'
145        if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
146            unittype='time'
147        elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc','months sin', 'years sinc']:
148            unittype='time'
149        elif string.lower(unit) in ['mbar', 'pa']:
150            unittype='pressure'
151        elif string.lower(unit) in ['m', 'km']:
152            unittype='height'
153        elif string.lower(unit) in ['degrees_north', 'degrees_south']:
154            unittype='latitude'
155        elif string.lower(unit) in ['degrees_east','degrees_west']:
156            unittype='longitude'
157        else:
158            unittype = unit
159        return unittype
160   
161    def determineCRS(self, axes=None, units=None, knownCRSAxes=None):
162        '''given any list of axis names and a list of units for these axes attempt to determine the CRS and return the CRSystem object and an axis order list.
163        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
164        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
165        Assumes the units are ordered to correspond with the axes
166        Alternatively, supply with a list of knownCRSAxes names (e.g. 'latitude', 'longitude', 'height') to get the right crs object       
167        '''     
168        axisorder=[]
169       
170        if knownCRSAxes is not None:
171            axes=knownCRSAxes
172        crs =None
173        #if len(axes)==2:
174            #crs=self.systems['ndg:crs:unknown2d']
175        #if len(axes)==3:
176            #crs=self.systems['ndg:crs:unknown3d']
177        #elif len(axes)==4:
178            #crs=self.systems['ndg:crs:unknown4d']
179        #else :
180            #crs=self.systems['ndg:crs:unknown']
181       
182        #this can be extended to accomodate more units       
183        if knownCRSAxes is not None:
184            crsMap=knownCRSAxes         
185        else:
186            crsMap=[]
187            for axis in axes:
188                unit=units[axes.index(axis)]
189                unittype=self.getUnitType(unit)           
190                crsMap.append(unittype)
191           
192       
193        #now try and match up crsMap with known crsystems.       
194        #axisorder stores the order of the supplied axes against the known crs.
195        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
196        #axisorder would be [2, 0, 1]         
197        for system in self.systems:
198            if len(self.systems[system].axes) != len(crsMap):
199                #wrong dimensionality, skip this crs
200                continue
201            match=0
202            axisorder={}
203            for ax in self.systems[system].axes:
204                if ax in crsMap:
205                    axisorder[ax]=crsMap.index(ax)
206                else:
207                    match = 1
208                    break
209            if match == 0: 
210                crs=self.systems[system]
211                break
212        if match==1:  #couldn't find a match - generate a pseudo crs           
213            for psystem in self.pseudosystems:
214                if len(self.pseudosystems[psystem].axes) == len(crsMap):
215                    axisorder={}
216                    indexes=[]
217                    for ax in crsMap:
218                        if ax in self.pseudosystems[psystem].axes:
219                            axisorder[ax]=crsMap.index(ax)
220                            indexes.append(crsMap.index(ax))
221                               
222                    if len(axisorder)== len(self.pseudosystems[psystem].axes) - 1:
223                        for ax in enumerate(crsMap):
224                            if ax[1] not in axisorder:
225                                for i in range(3):
226                                    if i not in indexes:                               
227                                        axisorder[axes[ax[0]]]=i
228                                        break         
229                                   
230                        axisLabels=''                 
231                        orderedNewAxes=[] 
232                        for ax in self.pseudosystems[psystem].axes:
233                            orderedNewAxes.append(ax)
234                        for key in axisorder:
235                            if key in orderedNewAxes:
236                                pass
237                            else:
238                                orderedNewAxes[orderedNewAxes.index('tba')] = key
239   
240                        for ax in orderedNewAxes:
241                            axisLabels=axisLabels + ax + ' '         
242                   
243                        crs=CRSystem('ndg:pseudoCRS', orderedNewAxes) #this should be ordered list.
244                        crs.axisLabels=axisLabels
245                        crs.timeAxis=self.pseudosystems[psystem].timeAxis
246                        crs.latAxis=self.pseudosystems[psystem].latAxis
247                        crs.latAxis=self.pseudosystems[psystem].lonAxis
248                        return crs, axisorder
249        return crs, axisorder
250       
251def main():
252    import cdms # just for testing
253    #sys.stdout = open("/home/dom/svn/trunk/csml/testfiles/gridseries/crsout.txt", "w")
254    cat=CRSCatalogue()
255 
256    #test every crs in this netcdf file
257    crsdictionary={}
258    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc')
259    for var in f.variables:
260        v=f.variables[var]
261        axes=v.getAxisIds()
262        units =[]
263        for axis in axes:
264            ax=f.getAxis(axis)
265            units.append(ax.attributes['units'])
266        crsdictionary[str(axes)]=units
267    n=0
268    #test each definition in dictionary
269       
270    for key in crsdictionary:
271        n=n+1
272        print '\n TEST %d'%n
273        print '*************************************'
274        #make list from string.
275        axs=[]
276        for item in key.split():
277            axs.append(item)
278        axs=eval(key)
279        uns=crsdictionary[key]
280        print 'AXES: %s'%axs
281        print 'UNITS: %s'%uns
282        print '----------------------------'
283        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
284        #print 'CRS: %s'%crs.srsName
285        print 'ORDER:  %s'%axisorder
286        print 'LABELS: %s'%crs.axisLabels
287        print 'CRS AXES: %s'%crs.axes
288        labels=''
289       
290       
291        for axis in axes:
292            labels=labels + axis + ' '
293       
294        #print '\n Pseudo CSML:'
295        #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
296        #print '<axisLabels>%s</axisLabels>'%labels
297        #for axis in axes:
298            #print '<GridOrdinate>'
299            #unit=units[axes.index(axis)]
300            #print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
301            #print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
302            #print '</GridOrdinate>'
303        #print '</GridSeriesDomain>'
304       
305        #print '**************************************'
306   #test getting CRS by name:
307    #print 'getting ndg:crs:xypt by name'
308    #crs=cat.getCRS('ndg:crs:xypt')
309    #print 'returns %s'%crs.srsName
310       
311   
312if __name__=="__main__":
313    main()
Note: See TracBrowser for help on using the repository browser.