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

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

handling missing units of measure

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