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

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

fixed bug that was caused by sorting axes list unintentionally

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