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

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

fix to CRS catalogue and to ops_gridseriesfeature

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