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

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

default calendar attribute added, repitition in fileList removed, gml:description fixed

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