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

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

adding pressure crs needed for BODC trajectories

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