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

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

added .twoD attribute to CRS sytems in catalogue - returns 2D element of nD systems if defined e.g. EPSG:4326

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