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

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

crs code working fully

Line 
1import string
2import sys         
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        self.pseudosystems={}
22       
23        # define lon lat pressure time CRS:       
24        crs=CRSystem(srsName='ndg:crs:xypt', axes =['longitude', 'latitude','pressure','time'])
25        crs.lonAxis=0
26        crs.latAxis=1
27        crs.timeAxis=3
28        self.systems['ndg:crs:xypt']=crs
29       
30        # define lon lat height time CRS:       
31        crs=CRSystem(srsName='ndg:crs:xyht', axes =['longitude', 'latitude','height','time'])
32        crs.lonAxis=0
33        crs.latAxis=1
34        crs.timeAxis=3
35        self.systems['ndg:crs:xyht']=crs
36       
37        # define lon lat time CRS:
38        crs=CRSystem(srsName='ndg:crs:xyt', axes =['longitude', 'latitude','time'])
39        self.systems['ndg:crs:xyt']=crs
40        crs.lonAxis=0
41        crs.latAxis=1
42        crs.timeAxis=2
43       
44        # define lon lat height CRS:
45        crs=CRSystem(srsName='ndg:crs:xyh', axes =['longitude', 'latitude','height'])
46        self.systems['ndg:crs:xyh']=crs
47        crs.lonAxis=0
48        crs.latAxis=1
49        crs.timeAxis=None
50       
51       
52        # define height time CRS:
53        crs=CRSystem(srsName='ndg:crs:ht', axes =['height', 'time'])
54        self.systems['ndg:crs:ht']=crs
55        crs.lonAxis=None
56        crs.latAxis=None
57        crs.timeAxis=1
58       
59        #define unknown 3D CRS:
60        crs=CRSystem(srsName='ndg:crs:undefined2d', axes=['tba','tba'])
61        self.pseudosystems['ndg:crs:undefined2d']=crs
62       
63        #define unknown 3D CRS:
64        crs=CRSystem(srsName='ndg:crs:undefinedll3d', axes=['longitude','latitude','tba'])
65        self.pseudosystems['ndg:crs:undefinedll3d']=crs
66        crs.lonAxis=0
67        crs.latAxis=1
68       
69        #define unknown 4D CRS:
70        crs=CRSystem(srsName='ndg:crs:undefinedllt4d', axes=['longitude','latitude','tba','time'])
71        self.pseudosystems['ndg:crs:undefinedllt4d']=crs
72        crs.lonAxis=0
73        crs.latAxis=1
74        crs.timeAxis=3
75       
76        #define unknown 1D CRS:
77        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
78        self.systems['ndg:crs:unknown']=crs
79                       
80        #define unknown 3D CRS:
81        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
82        self.systems['ndg:crs:unknown2d']=crs
83       
84        #define unknown 3D CRS:
85        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
86        self.systems['ndg:crs:unknown3d']=crs
87       
88        #define unknown 4D CRS:
89        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
90        self.systems['ndg:crs:unknown4d']=crs
91       
92           
93    def getCRS(self, srsName, labels=None):
94        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
95        #labels are only required if it is a ndg:pseudoCRS
96        try:
97            return self.systems[srsName]
98        except:
99            if srsName=='ndg:pseudoCRS':
100                labellist=[]
101                for label in labels.split():
102                    labellist.append(label)
103            #need to find the right pseudosystem.
104            #at the moment a simple test on number of dimensions works, but need to make this more generic.
105            for psystem in self.pseudosystems:
106                if len(self.pseudosystems[psystem].axes) == len (labellist):
107                    return self.pseudosystems[psystem]
108            else:
109                return self.systems['ndg:crs:unknown']
110   
111    def getUnitType(self, unit):
112        #unittype='unknown'
113        if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
114            unittype='time'
115        elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc','months sin', 'years sinc']:
116            unittype='time'
117        elif string.lower(unit) in ['mbar', 'pa']:
118            unittype='pressure'
119        elif string.lower(unit) in ['m', 'km']:
120            unittype='height'
121        elif string.lower(unit) in ['degrees_north', 'degrees_south']:
122            unittype='latitude'
123        elif string.lower(unit) in ['degrees_east','degrees_west']:
124            unittype='longitude'
125        else:
126            unittype = unit
127        return unittype
128   
129    def determineCRS(self, axes=None, units=None, knownCRSAxes=None):
130        '''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.
131        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
132        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
133        Assumes the units are ordered to correspond with the axes
134        Alternatively, supply with a list of knownCRSAxes names (e.g. 'latitude', 'longitude', 'height') to get the right crs object       
135        '''     
136        axisorder=[]
137
138        if knownCRSAxes is not None:
139            axes=knownCRSAxes
140        crs =None
141        #if len(axes)==2:
142            #crs=self.systems['ndg:crs:unknown2d']
143        #if len(axes)==3:
144            #crs=self.systems['ndg:crs:unknown3d']
145        #elif len(axes)==4:
146            #crs=self.systems['ndg:crs:unknown4d']
147        #else :
148            #crs=self.systems['ndg:crs:unknown']
149       
150        #this can be extended to accomodate more units       
151        if knownCRSAxes is not None:
152            crsMap=knownCRSAxes         
153        else:
154            crsMap=[]
155            for axis in axes:
156                unit=units[axes.index(axis)]
157                unittype=self.getUnitType(unit)           
158                crsMap.append(unittype)
159           
160       
161        #now try and match up crsMap with known crsystems.       
162        #axisorder stores the order of the supplied axes against the known crs.
163        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
164        #axisorder would be [2, 0, 1]         
165        for system in self.systems:
166            if len(self.systems[system].axes) != len(crsMap):
167                #wrong dimensionality, skip this crs
168                continue
169            match=0
170            axisorder={}
171            for ax in self.systems[system].axes:
172                if ax in crsMap:
173                    axisorder[ax]=crsMap.index(ax)
174                else:
175                    match = 1
176                    break
177            if match == 0: 
178                crs=self.systems[system]
179                break
180        if match==1:  #couldn't find a match - generate a pseudo crs   
181            for psystem in self.pseudosystems:
182                if len(self.pseudosystems[psystem].axes) != len(crsMap):
183                    #wrong dimensionality, skip this crs
184                    continue
185                axisorder={}
186                indexes=[]
187                for ax in crsMap:
188                    if ax in self.pseudosystems[psystem].axes:
189                        axisorder[ax]=crsMap.index(ax)
190                        indexes.append(crsMap.index(ax))
191                             
192                if len(axisorder)== len(self.pseudosystems[psystem].axes) - 1:
193                    for ax in enumerate(crsMap):
194                        if ax[1] not in axisorder:
195                            for i in range(3):
196                                if i not in indexes:                               
197                                    axisorder[axes[ax[0]]]=i
198                                    break         
199                               
200                    axisLabels=''                 
201                    orderedNewAxes=[] 
202                    for ax in self.pseudosystems[psystem].axes:
203                        orderedNewAxes.append(ax)
204                    for key in axisorder:
205                        if key in orderedNewAxes:
206                            pass
207                        else:
208                            orderedNewAxes[orderedNewAxes.index('tba')] = key
209
210                    for ax in orderedNewAxes:
211                        axisLabels=axisLabels + ax + ' '         
212                 
213                    crs=CRSystem('ndg:pseudoCRS', orderedNewAxes) #this should be ordered list.
214                    crs.axisLabels=axisLabels
215                    crs.timeAxis=self.pseudosystems[psystem].timeAxis
216                    crs.latAxis=self.pseudosystems[psystem].latAxis
217                    crs.latAxis=self.pseudosystems[psystem].lonAxis
218        return crs, axisorder
219       
220def main():
221    import cdms # just for testing
222    #sys.stdout = open("/home/dom/svn/trunk/csml/testfiles/gridseries/crsout.txt", "w")
223    cat=CRSCatalogue()
224 
225    #test every crs in this netcdf file
226    crsdictionary={}
227    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc')
228    for var in f.variables:
229        v=f.variables[var]
230        axes=v.getAxisIds()
231        units =[]
232        for axis in axes:
233            ax=f.getAxis(axis)
234            units.append(ax.attributes['units'])
235        crsdictionary[str(axes)]=units
236    n=0
237    #test each definition in dictionary
238       
239    for key in crsdictionary:
240        n=n+1
241        print '\n TEST %d'%n
242        print '*************************************'
243        #make list from string.
244        axs=[]
245        for item in key.split():
246            axs.append(item)
247        axs=eval(key)
248        uns=crsdictionary[key]
249        print 'AXES: %s'%axs
250        print 'UNITS: %s'%uns
251        print '----------------------------'
252        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
253        #print 'CRS: %s'%crs.srsName
254        print 'ORDER:  %s'%axisorder
255        print 'LABELS: %s'%crs.axisLabels
256        print 'CRS AXES: %s'%crs.axes
257        labels=''
258       
259       
260        for axis in axes:
261            labels=labels + axis + ' '
262       
263        #print '\n Pseudo CSML:'
264        #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
265        #print '<axisLabels>%s</axisLabels>'%labels
266        #for axis in axes:
267            #print '<GridOrdinate>'
268            #unit=units[axes.index(axis)]
269            #print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
270            #print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
271            #print '</GridOrdinate>'
272        #print '</GridSeriesDomain>'
273       
274        #print '**************************************'
275   #test getting CRS by name:
276    #print 'getting ndg:crs:xypt by name'
277    #crs=cat.getCRS('ndg:crs:xypt')
278    #print 'returns %s'%crs.srsName
279       
280   
281if __name__=="__main__":
282    main()
Note: See TracBrowser for help on using the repository browser.