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

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

removing obsolete axisList class

Line 
1import string
2import cdms # just while testing
3
4class CRSystem(object):
5    def __init__(self, srsName, axes):
6        self.srsName=srsName
7        self.axes=axes
8        self.srsDimension=len(self.axes)
9        self.axisLabels=''
10        for axis in self.axes:
11            self.axisLabels=self.axisLabels + axis + ' '
12        self.timeAxis=None
13        self.lonAxis=None
14        self.latAxis=None
15
16class CRSCatalogue(object):
17    def __init__(self):
18        # dictionary to hold known CRSystems:
19        self.systems={}
20       
21        # define lon lat pressure time CRS:       
22        crs=CRSystem(srsName='ndg:crs:xypt', axes =['Lon', 'Lat','Pressure','Time'])
23        crs.lonAxis=0
24        crs.latAxis=1
25        crs.timeAxis=3
26        self.systems['ndg:crs:xypt']=crs
27       
28        # define lon lat height time CRS:       
29        crs=CRSystem(srsName='ndg:crs:xyht', axes =['Lon', 'Lat','Height','Time'])
30        crs.lonAxis=0
31        crs.latAxis=1
32        crs.timeAxis=3
33        self.systems['ndg:crs:xyht']=crs
34       
35        # define lon lat time CRS:
36        crs=CRSystem(srsName='ndg:crs:xyt', axes =['Lon', 'Lat','Time'])
37        self.systems['ndg:crs:xyt']=crs
38       
39        #define unknown 1D CRS:
40        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
41        self.systems['ndg:crs:unknown']=crs
42       
43        #define unknown 3D CRS:
44        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
45        self.systems['ndg:crs:unknown2d']=crs
46       
47        #define unknown 3D CRS:
48        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
49        self.systems['ndg:crs:unknown3d']=crs
50       
51        #define unknown 4D CRS:
52        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
53        self.systems['ndg:crs:unknown4d']=crs
54       
55           
56    def getCRS(self, srsName):
57        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
58        try:
59            return self.systems[srsName]
60        except:
61            return self.systems['ndg:crs:unknown']
62   
63    def determineCRS(self, axes, units):
64        '''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.
65        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
66        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
67        Assumes the units are ordered to correspond with the axes!
68        '''
69        axisorder=[]
70        if len(axes)==2:
71            crs=self.systems['ndg:crs:unknown2d']
72        if len(axes)==3:
73            crs=self.systems['ndg:crs:unknown3d']
74        elif len(axes)==4:
75            crs=self.systems['ndg:crs:unknown4d']
76        else :
77            crs=self.systems['ndg:crs:unknown']
78        crsMap=[]
79        #this can be extended to accomodate more units
80        for axis in axes:
81            unit=units[axes.index(axis)]
82            unittype='unknown'
83            if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
84                unittype='Time'
85            elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc']:
86                unittype='Time'
87            elif string.lower(unit) in ['mbar', 'pa','level']:
88                unittype='Pressure'
89            elif string.lower(unit) in ['m', 'km']:
90                unittype='Height'
91            elif string.lower(unit) in ['degrees_north']:
92                unittype='Lat'
93            elif string.lower(unit) in ['degrees_east']:
94                unittype='Lon'
95            crsMap.append(unittype)
96       
97        #now try and match up crsMap with known crsystems.       
98        #axisorder stores the order of the supplied axes against the known crs.
99        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
100        #axisorder would be [2, 0, 1] 
101       
102        for system in self.systems:
103            if len(self.systems[system].axes) != len(crsMap):
104                #wrong dimensionality, skip this crs
105                continue
106            match=0
107            axisorder=[]
108            for ax in self.systems[system].axes:
109                if ax in crsMap:
110                    axisorder.append(crsMap.index(ax))
111                else:
112                    match = 1
113                    break
114            if match == 0: 
115                crs=self.systems[system]
116                break
117                       
118        return crs, axisorder
119       
120def main():
121    cat=CRSCatalogue()
122 
123
124    #test every crs in this netcdf file
125    crsdictionary={}
126    f=cdms.open('/home/dom/Desktop/SVN/trunk/csml/testfiles/gridseries/xaaqda@pxs19c1.nc')
127    for var in f.variables:
128        v=f.variables[var]
129        axes=v.getAxisIds()
130        units =[]
131        for axis in axes:
132            ax=f.getAxis(axis)
133            units.append(ax.attributes['units'])
134        crsdictionary[str(axes)]=units
135    n=2
136    #test each definition in dictionary
137    for key in crsdictionary:
138        n=n+1
139        print '\n TEST %d'%n
140        print '*************************************'
141       #make list from string.
142        axs=[]
143        for item in key.split():
144            axs.append(item)
145        uns=crsdictionary[key]
146        print 'AXES: %s'%axs
147        print 'UNITS: %s'%uns
148        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
149        print 'CRS: %s'%crs.srsName
150        print 'ORDER:  %s'%axisorder
151        #print 'axis order =%s'%axisorder
152        labels=''
153       
154        print axes
155        # ['t', 'ht', 'latitude_1', 'longitude_1']
156        print crs.axes
157         # ['Lon', 'Lat', 'Pressure', 'Time']
158        print axisorder
159        # [3, 2, 1, 0]
160       
161        for axis in axes:
162            labels=labels + axis + ' '
163       
164        print '\n Pseudo CSML:'
165        print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
166        print '<axisLabels>%s</axisLabels>'%labels
167        for axis in axes:
168            print '<GridOrdinate>'
169 
170            unit=units[axes.index(axis)]
171            print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
172            print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
173            print '</GridOrdinate>'
174        print '</GridSeriesDomain>'
175       
176        print '**************************************'
177   #test getting CRS by name:
178    #print 'getting ndg:crs:xypt by name'
179    #crs=cat.getCRS('ndg:crs:xypt')
180    #print 'returns %s'%crs.srsName
181       
182   
183if __name__=="__main__":
184    main()
Note: See TracBrowser for help on using the repository browser.