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

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

added new crs definition to catalogue and also exception raised when suitable definition is not found

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