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

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

rewriting of code that determines CRS system

Line 
1import string
2import sys     
3import csml.csmllibs.csmlextra
4from copy import copy
5
6class CRSystem(object):
7    def __init__(self, srsName, axes):
8        self.srsName=srsName
9        self.axes=axes
10        self.srsDimension=len(self.axes)
11        self.axisLabels=''
12        for axis in self.axes:
13            self.axisLabels=self.axisLabels + axis + ' '
14        self.axisLabels=self.axisLabels[:-1]
15        self.timeAxis=None
16        self.lonAxis=None
17        self.latAxis=None
18        self.units=None
19   
20    def getAxisLabels(self):
21        labelList=[]
22        for item in self.axisLabels.split():
23            labelList.append(item)
24        return labelList
25
26class CRSCatalogue(object):
27    def __init__(self):
28        # dictionary to hold known CRSystems:
29        self.systems={}
30        self.pseudosystems={}
31       
32        # define WGS84 CRS:        #time?
33        crs=CRSystem(srsName='WGS84', axes =['time', 'latitude','longitude'])
34        crs.lonAxis=2
35        crs.latAxis=1
36        crs.timeAxis=0
37        self.systems['WGS84']=crs
38       
39        # define WGS84 CRS:        #no time?
40        #crs=CRSystem(srsName='WGS84', axes =['longitude','latitude'])
41        #crs.lonAxis=0
42        #crs.latAxis= 1
43        ##crs.timeAxis=0
44        #self.systems['WGS84']=crs
45       
46       
47         #time only crs
48        crs=CRSystem(srsName='ndg:crs:t', axes =['time'])
49        crs.timeAxis=0
50        crs.lonAxis=None
51        crs.latAxis=None
52        self.systems['ndg:crs:t']=crs
53       
54        #pressure only crs
55        crs=CRSystem(srsName='ndg:crs:p', axes =['pressure'])
56        crs.timeAxis=0
57        crs.lonAxis=None
58        crs.latAxis=None
59        self.systems['ndg:crs:p']=crs
60       
61       
62        # define lon lat pressure time CRS:       
63        crs=CRSystem(srsName='ndg:crs:lonlatpt', axes =['longitude', 'latitude','pressure','time'])
64        crs.lonAxis=0
65        crs.latAxis=1
66        crs.timeAxis=3
67        self.systems['ndg:crs:lonlatpt']=crs
68       
69        # define lon lat height time CRS:       
70        crs=CRSystem(srsName='ndg:crs:lonlatht', axes =['longitude', 'latitude','height','time'])
71        crs.lonAxis=0
72        crs.latAxis=1
73        crs.timeAxis=3
74        self.systems['ndg:crs:lonlatht']=crs
75       
76        # define lon lat time CRS:
77        crs=CRSystem(srsName='ndg:crs:lonlatt', axes =['longitude', 'latitude','time'])
78        self.systems['ndg:crs:lonlatt']=crs
79        crs.lonAxis=0
80        crs.latAxis=1
81        crs.timeAxis=2
82       
83        # define lon lat height CRS:
84        crs=CRSystem(srsName='ndg:crs:lonlath', axes =['longitude', 'latitude','height'])
85        self.systems['ndg:crs:lonlath']=crs
86        crs.lonAxis=0
87        crs.latAxis=1
88        crs.timeAxis=None
89       
90       
91        # define height time CRS:
92        crs=CRSystem(srsName='ndg:crs:ht', axes =['height', 'time'])
93        self.systems['ndg:crs:ht']=crs
94        crs.lonAxis=None
95        crs.latAxis=None
96        crs.timeAxis=1
97       
98        # define british national grid plus time CRS:
99        crs=CRSystem(srsName='ndg:crs:bngt', axes =['eastings','northings','time'])
100        self.systems['ndg:crs:bngt']=crs
101        crs.lonAxis=0
102        crs.latAxis=1
103        crs.timeAxis=2
104       
105       
106        #define rotated_latitude_longitude CRS (test)
107        #crs=CRSystem(srsName='ndg:crs:rlonrlatsurft', axes =['grid_longitude', 'grid_latitude','surface', 'time'])
108        #self.systems['ndg:crs:rlonrlatsurft']=crs
109        #crs.lonAxis=0
110        #crs.latAxis=1
111        #crs.timeAxis=3
112       
113       
114        ##define rotated_latitude_longitude_undefined surface CRS (test)
115        #crs=CRSystem(srsName='ndg:crs:undefined_rlonrlatt', axes =['grid_longitude', 'grid_latitude','tba', 'time'])
116        #self.pseudosystems['ndg:crs:undefined_rlonrlatt']=crs
117        #crs.lonAxis=0
118        #crs.latAxis=1
119        #crs.timeAxis=3
120       
121        #define unknown 2D CRS:
122        crs=CRSystem(srsName='ndg:crs:undefined_xt', axes=['tba','time'])
123        self.pseudosystems['ndg:crs:undefined_xt']=crs
124        crs.timeAxis=1
125        crs.lonAxis=None
126        crs.latAxis=None
127       
128        #define unknown 3D CRS:
129        crs=CRSystem(srsName='ndg:crs:undefined_lonlatz', axes=['longitude','latitude','tba'])
130        self.pseudosystems['ndg:crs:undefined_lonlatz']=crs
131        crs.lonAxis=0
132        crs.latAxis=1
133        crs.timeAxis=None
134       
135        #define unknown 4D CRS:
136        crs=CRSystem(srsName='ndg:crs:undefined_lonlatzt', axes=['longitude','latitude','tba','time'])
137        self.pseudosystems['ndg:crs:undefined_lonlatzt']=crs
138        crs.lonAxis=0
139        crs.latAxis=1
140        crs.timeAxis=3
141       
142        #define unknown 5D CRS (test):
143        crs=CRSystem(srsName='ndg:crs:undefined5d_lonlatzt', axes=['probAx','time','tba', 'latitude','longitude'])
144        self.pseudosystems['ndg:crs:undefined5d_lonlatzt']=crs
145        crs.lonAxis=4
146        crs.latAxis=3
147        crs.timeAxis=1
148       
149        #define unknown 1D CRS:
150        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
151        self.systems['ndg:crs:unknown']=crs
152                       
153        #define unknown 3D CRS:
154        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
155        self.systems['ndg:crs:unknown2d']=crs
156       
157        #define unknown 3D CRS:
158        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
159        self.systems['ndg:crs:unknown3d']=crs
160       
161        #define unknown 4D CRS:
162        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
163        self.systems['ndg:crs:unknown4d']=crs
164                   
165    def getCRS(self, srsName):
166        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
167        try:
168            return self.systems[srsName]
169        except:
170            return self.pseudosystems[srsName]
171   
172   
173    def getUnitType(self, unit):
174        #unittype='unknown'
175        if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
176            unittype='time'
177        elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc','months sin', 'years sinc']:
178            unittype='time'
179        elif string.lower(unit) in ['mbar', 'pa']:
180            unittype='pressure'
181        elif string.lower(unit) in ['m', 'km']:
182            unittype='height'
183        elif string.lower(unit) in ['degrees_north', 'degrees_south']:
184            unittype='latitude'
185        elif string.lower(unit) in ['degrees_east','degrees_west']:
186            unittype='longitude'
187        elif string.lower(unit) in ['-']:
188            unittype='probAx'
189           
190        else:
191            unittype = unit
192        return unittype
193
194    def __generatePseudoCRS(self,crsMap, axes, units):
195        ''' if a suitable crs cannot be found a pseudo crs is created'''
196        for psystem in self.pseudosystems:
197            if len(self.pseudosystems[psystem].axes) == len(crsMap):
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=CRSystem(psystem, orderedNewAxes) #this should be ordered list.
228                    crs.axisLabels=axisLabels
229                    crs.timeAxis=self.pseudosystems[psystem].timeAxis
230                    crs.latAxis=self.pseudosystems[psystem].latAxis
231                    crs.latAxis=self.pseudosystems[psystem].lonAxis
232                    orderedunits=[]
233                    for a in crs.axes:
234                        orderedunits.append(units[axisorder[a]])     
235                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits)                         
236        return crs, axisorder
237
238   
239    def _compareRealAxisNames(self, axes):
240        ''' will find complete exact matches only'''
241        crs = None
242        axisorder={}
243        unsortedaxes=copy(axes)
244        for system in self.systems:
245            crsaxes=self.systems[system].axes
246            unsortedcrsaxes=copy(crsaxes)
247            crsaxes.sort()
248            axes.sort()
249            if crsaxes==axes:
250                #found exact match  - stop checking
251                crs=self.systems[system]
252                axisorder=self._getAxisOrder(unsortedaxes,unsortedcrsaxes)
253                break
254        return crs, axisorder
255       
256    def _getAxisOrder(self, axes, crsaxes):
257        #given a list of axis names and a list of comparable crsaxes returns an axis order list e.g. [3,2,1,0] where '3' means the first item in axes corresponds to the item at index 3 in the crsaxes list.
258        axisorder={}
259        for ax in axes:
260                if ax in crsaxes:
261                    axisorder[ax]=crsaxes.index(ax)
262        return axisorder
263   
264    def _compareNamesAndUnits(self, axes, units, crsMap, stdNames=None):
265        #build a map of axes and units - use std names to help!
266        crs = None
267        axisorder={}
268        if not crsMap:
269            for axis in axes:
270                unittype = None
271                if stdNames is not None:
272                    if stdNames[axes.index(axis)]:
273                        unittype=stdNames[axes.index(axis)]
274                else:
275                    unit=units[axes.index(axis)]
276                    unittype=self.getUnitType(unit)         
277                crsMap.append(unittype) 
278           
279        match=0
280        #now try and match up crsMap with known crsystems.       
281        #axisorder stores the order of the supplied axes against the known crs.
282        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
283        #axisorder would be [2, 0, 1]         
284        for system in self.systems:
285            if len(self.systems[system].axes) != len(crsMap):
286                #wrong dimensionality, skip this crs
287                continue
288            match=1
289            axisorder={}
290            for ax in self.systems[system].axes:
291                if ax in crsMap:
292                    axisorder[ax]=crsMap.index(ax)
293                else:
294                    match = 0
295                    break
296            if match == 1: 
297                crs=self.systems[system]             
298                break   
299        return crs, axisorder, crsMap
300   
301    def _addUnitsToCRSDefinition(self, crs, units, axisorder):
302        '''if no units are defined in the crs, set the units to be those of the original data'''
303        if crs.units is None:
304            if units is not None:
305                orderedunits=[]
306                for a in crs.axes:
307                    orderedunits.append(units[axisorder[a]])
308                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits)
309        return crs
310               
311    def determineCRS(self, axes=None, units=None, knownCRSAxes=None,stdNames=None):
312        '''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.
313        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
314        returns self.systems['ndg:crs:xypt'], {'latitude': 2, 'pressure': 1, 'longitude': 3, 'time': 0}
315        Assumes the units are ordered to correspond with the axes
316        ''' 
317        axisorder={}
318        crs =None
319        #perform various tests to try and determine which crs this is.
320       
321        #test: compare real names of axes with names of axes in coordinate reference systems
322        #looks for exact matches only.
323        if axes is not None:
324            crs, axisorder = self._compareRealAxisNames(copy(axes))
325       
326        #set up crsMap - may be empty
327        crsMap=[]
328        if knownCRSAxes is not None:
329            for item in knownCRSAxes:
330                if item is not None:
331                    crsMap=knownCRSAxes
332                    break 
333       
334       
335        if crs is None:
336            #try and match based on units, names and stdNames
337            crs,axisorder, crsMap=self._compareNamesAndUnits(axes, units, crsMap, stdNames=None)
338        if crs is None:
339            #try to generate a pseudo crs
340            crs, axisorder=self.__generatePseudoCRS(crsMap, axes, units)
341           
342        #Finally, if no units are defined in the crs, set the units to be those of the original data
343        if crs is not None:
344            crs=self._addUnitsToCRSDefinition(crs, units, axisorder)
345        return crs, axisorder
346     
347    def _identifyAxisByStdName(self,  stdName):
348        '''given a CF standard name return the axis name'''
349        understoodNames={'grid_longitude': 'grid_longitude', 'grid_latitude':'grid_latitude'}
350        if stdName in understoodNames.keys():
351            axisname= understoodNames[stdName]
352        else:
353            axisname =None
354        return axisname
355       
356def main():
357    import cdms # just for testing
358    #sys.stdout = open("/home/dom/svn/trunk/csml/testfiles/gridseries/crsout.txt", "w")
359    cat=CRSCatalogue()
360 
361    #test every crs in this netcdf file
362    crsdictionary={}
363    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc')
364    #f=cdms.open('/home/dom/EN/air_frost_1961-1990.nc')
365    for var in f.variables:
366        v=f.variables[var]
367        axes=v.getAxisIds()
368        units =[]
369        for axis in axes:
370            ax=f.getAxis(axis)
371            units.append(ax.attributes['units'])
372        crsdictionary[str(axes)]=units
373    n=0
374    #test each definition in dictionary
375       
376    for key in crsdictionary:
377        n=n+1
378        print '\n TEST %d'%n
379        print '*************************************'
380        #make list from string.
381        axs=[]
382        for item in key.split():
383            axs.append(item)
384        axs=eval(key)
385        uns=crsdictionary[key]
386        print 'AXES: %s'%axs
387        print 'UNITS: %s'%uns
388        print '----------------------------'
389        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
390        print crs, axisorder
391        if crs is not None:
392            print 'CRS found: %s'%crs.srsName
393            print 'ORDER:  %s'%axisorder
394            print 'LABELS: %s'%crs.axisLabels
395            print 'CRS AXES: %s'%crs.axes
396        else:
397            print 'CRS not found'
398        labels=''
399       
400       
401        #for axis in axes:
402            #labels=labels + axis + ' '
403       
404        ##print '\n Pseudo CSML:'
405        ##print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
406        ##print '<axisLabels>%s</axisLabels>'%labels
407        ##for axis in axes:
408            ##print '<GridOrdinate>'
409            ##unit=units[axes.index(axis)]
410            ##print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
411            ##print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
412            ##print '</GridOrdinate>'
413        ##print '</GridSeriesDomain>'
414       
415        ##print '**************************************'
416   ##test getting CRS by name:
417    #print 'getting ndg:crs:xypt by name'
418    #crs=cat.getCRS('ndg:crs:xypt')
419    #print 'returns %s'%crs.srsName
420       
421   
422if __name__=="__main__":
423    main()
Note: See TracBrowser for help on using the repository browser.