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

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

some CRS issues resolved

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