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

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

added ability to filter out low dimension variables in scanner. Added some test crs definitions for testing rotated grids. Added axis identification by standard name

Line 
1import string
2import sys     
3import csml.csmllibs.csmlextra
4
5class CRSystem(object):
6    def __init__(self, srsName, axes):
7        self.srsName=srsName
8        self.axes=axes
9        self.srsDimension=len(self.axes)
10        self.axisLabels=''
11        for axis in self.axes:
12            self.axisLabels=self.axisLabels + axis + ' '
13        self.axisLabels=self.axisLabels[:-1]
14        self.timeAxis=None
15        self.lonAxis=None
16        self.latAxis=None
17        self.units=None
18   
19    def getAxisLabels(self):
20        labelList=[]
21        for item in self.axisLabels.split():
22            labelList.append(item)
23        return labelList
24
25class CRSCatalogue(object):
26    def __init__(self):
27        # dictionary to hold known CRSystems:
28        self.systems={}
29        self.pseudosystems={}
30       
31        # define WGS84 CRS:        #time?
32        crs=CRSystem(srsName='WGS84', axes =['time', 'latitude','longitude'])
33        crs.lonAxis=2
34        crs.latAxis=1
35        crs.timeAxis=0
36        self.systems['WGS84']=crs
37       
38        # define WGS84 CRS:        #no time?
39        #crs=CRSystem(srsName='WGS84', axes =['longitude','latitude'])
40        #crs.lonAxis=0
41        #crs.latAxis= 1
42        ##crs.timeAxis=0
43        #self.systems['WGS84']=crs
44       
45       
46         #time only crs
47        crs=CRSystem(srsName='ndg:crs:t', axes =['time'])
48        crs.timeAxis=0
49        crs.lonAxis=None
50        crs.latAxis=None
51        self.systems['ndg:crs:t']=crs
52       
53        #pressure only crs
54        crs=CRSystem(srsName='ndg:crs:p', axes =['pressure'])
55        crs.timeAxis=0
56        crs.lonAxis=None
57        crs.latAxis=None
58        self.systems['ndg:crs:p']=crs
59       
60       
61        # define lon lat pressure time CRS:       
62        crs=CRSystem(srsName='ndg:crs:lonlatpt', axes =['longitude', 'latitude','pressure','time'])
63        crs.lonAxis=0
64        crs.latAxis=1
65        crs.timeAxis=3
66        self.systems['ndg:crs:lonlatpt']=crs
67       
68        # define lon lat height time CRS:       
69        crs=CRSystem(srsName='ndg:crs:lonlatht', axes =['longitude', 'latitude','height','time'])
70        crs.lonAxis=0
71        crs.latAxis=1
72        crs.timeAxis=3
73        self.systems['ndg:crs:lonlatht']=crs
74       
75        # define lon lat time CRS:
76        crs=CRSystem(srsName='ndg:crs:lonlatt', axes =['longitude', 'latitude','time'])
77        self.systems['ndg:crs:lonlatt']=crs
78        crs.lonAxis=0
79        crs.latAxis=1
80        crs.timeAxis=2
81       
82        # define lon lat height CRS:
83        crs=CRSystem(srsName='ndg:crs:lonlath', axes =['longitude', 'latitude','height'])
84        self.systems['ndg:crs:lonlath']=crs
85        crs.lonAxis=0
86        crs.latAxis=1
87        crs.timeAxis=None
88       
89       
90        # define height time CRS:
91        crs=CRSystem(srsName='ndg:crs:ht', axes =['height', 'time'])
92        self.systems['ndg:crs:ht']=crs
93        crs.lonAxis=None
94        crs.latAxis=None
95        crs.timeAxis=1
96       
97        #define rotated_latitude_longitude CRS (test)
98        crs=CRSystem(srsName='ndg:crs:rlonrlatsurft', axes =['grid_longitude', 'grid_latitude','surface', 'time'])
99        self.systems['ndg:crs:rlonrlatsurft']=crs
100        crs.lonAxis=0
101        crs.latAxis=1
102        crs.timeAxis=3
103       
104       
105        #define rotated_latitude_longitude_undefined surface CRS (test)
106        crs=CRSystem(srsName='ndg:crs:undefined_rlonrlatt', axes =['grid_longitude', 'grid_latitude','tba', 'time'])
107        self.pseudosystems['ndg:crs:undefined_rlonrlatt']=crs
108        crs.lonAxis=0
109        crs.latAxis=1
110        crs.timeAxis=3
111       
112        #define unknown 2D CRS:
113        crs=CRSystem(srsName='ndg:crs:undefined_xt', axes=['tba','time'])
114        self.pseudosystems['ndg:crs:undefined_xt']=crs
115        crs.timeAxis=1
116        crs.lonAxis=None
117        crs.latAxis=None
118       
119        #define unknown 3D CRS:
120        crs=CRSystem(srsName='ndg:crs:undefined_lonlatz', axes=['longitude','latitude','tba'])
121        self.pseudosystems['ndg:crs:undefined_lonlatz']=crs
122        crs.lonAxis=0
123        crs.latAxis=1
124        crs.timeAxis=None
125       
126        #define unknown 4D CRS:
127        crs=CRSystem(srsName='ndg:crs:undefined_lonlatzt', axes=['longitude','latitude','tba','time'])
128        self.pseudosystems['ndg:crs:undefined_lonlatzt']=crs
129        crs.lonAxis=0
130        crs.latAxis=1
131        crs.timeAxis=3
132       
133        #define unknown 5D CRS (test):
134        crs=CRSystem(srsName='ndg:crs:undefined5d_lonlatzt', axes=['probAx','time','tba', 'latitude','longitude'])
135        self.pseudosystems['ndg:crs:undefined5d_lonlatzt']=crs
136        crs.lonAxis=4
137        crs.latAxis=3
138        crs.timeAxis=1
139       
140        #define unknown 1D CRS:
141        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
142        self.systems['ndg:crs:unknown']=crs
143                       
144        #define unknown 3D CRS:
145        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
146        self.systems['ndg:crs:unknown2d']=crs
147       
148        #define unknown 3D CRS:
149        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
150        self.systems['ndg:crs:unknown3d']=crs
151       
152        #define unknown 4D CRS:
153        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
154        self.systems['ndg:crs:unknown4d']=crs
155                   
156    def getCRS(self, srsName):
157        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
158        try:
159            return self.systems[srsName]
160        except:
161            return self.pseudosystems[srsName]
162   
163   
164    def getUnitType(self, unit):
165        #unittype='unknown'
166        if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
167            unittype='time'
168        elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc','months sin', 'years sinc']:
169            unittype='time'
170        elif string.lower(unit) in ['mbar', 'pa']:
171            unittype='pressure'
172        elif string.lower(unit) in ['m', 'km']:
173            unittype='height'
174        elif string.lower(unit) in ['degrees_north', 'degrees_south']:
175            unittype='latitude'
176        elif string.lower(unit) in ['degrees_east','degrees_west']:
177            unittype='longitude'
178        elif string.lower(unit) in ['-']:
179            unittype='probAx'
180           
181        else:
182            unittype = unit
183        return unittype
184
185    def __generatePseudoCRS(self,crsMap, axes, units):
186        ''' if a suitable crs cannot be found a pseudo crs is created'''
187        for psystem in self.pseudosystems:
188            if len(self.pseudosystems[psystem].axes) == len(crsMap):
189                axisorder={}
190                indexes=[]
191                for ax in crsMap:
192                    if ax in self.pseudosystems[psystem].axes:
193                        axisorder[ax]=crsMap.index(ax)
194                        indexes.append(crsMap.index(ax))
195                           
196                if len(axisorder)== len(self.pseudosystems[psystem].axes) - 1:
197                    for ax in enumerate(crsMap):
198                        if ax[1] not in axisorder:
199                            for i in range(3):
200                                if i not in indexes:                               
201                                    axisorder[axes[ax[0]]]=i
202                                    break         
203                               
204                    axisLabels=''                 
205                    orderedNewAxes=[] 
206                    for ax in self.pseudosystems[psystem].axes:
207                        orderedNewAxes.append(ax)
208                    for key in axisorder:
209                        if key in orderedNewAxes:
210                            pass
211                        else:
212                            orderedNewAxes[orderedNewAxes.index('tba')] = key
213
214                    for ax in orderedNewAxes:
215                        axisLabels=axisLabels + ax + ' '         
216               
217                    #crs=CRSystem('ndg:pseudoCRS', orderedNewAxes) #this should be ordered list.
218                    crs=CRSystem(psystem, orderedNewAxes) #this should be ordered list.
219                    crs.axisLabels=axisLabels
220                    crs.timeAxis=self.pseudosystems[psystem].timeAxis
221                    crs.latAxis=self.pseudosystems[psystem].latAxis
222                    crs.latAxis=self.pseudosystems[psystem].lonAxis
223                    orderedunits=[]
224                    for a in crs.axes:
225                        orderedunits.append(units[axisorder[a]])     
226                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits)                         
227        return crs, axisorder
228
229    def determineCRS(self, axes=None, units=None, stdNames=None):
230        '''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.
231        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
232        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
233        Assumes the units are ordered to correspond with the axes
234        '''     
235        axisorder=[]
236        crs =None
237       
238        #build a map of axes and units - use std names to help!
239        crsMap=[]
240        for axis in axes:
241            if stdNames[axes.index(axis)]:
242                unittype=stdNames[axes.index(axis)]
243            else:
244                unit=units[axes.index(axis)]
245                unittype=self.getUnitType(unit)           
246            crsMap.append(unittype)
247             
248        match=0
249        #now try and match up crsMap with known crsystems.       
250        #axisorder stores the order of the supplied axes against the known crs.
251        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
252        #axisorder would be [2, 0, 1]         
253        for system in self.systems:
254            if len(self.systems[system].axes) != len(crsMap):
255                #wrong dimensionality, skip this crs
256                continue
257            match=1
258            axisorder={}
259            for ax in self.systems[system].axes:
260                if ax in crsMap:
261                    axisorder[ax]=crsMap.index(ax)
262                else:
263                    match = 0
264                    break
265            if match == 1: 
266                crs=self.systems[system]             
267                break
268       
269        if match==0:  #if couldn't find a match - generate a pseudo crs           
270            crs, axisorder=self.__generatePseudoCRS(crsMap, axes, units)
271       
272        if crs.units is None:
273            #if no units are defined, set the units to be those of the original data
274            if units is not None:
275                orderedunits=[]
276                for a in crs.axes:
277                    orderedunits.append(units[axisorder[a]])
278                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits)
279        return crs, axisorder
280   
281    def _identifyAxisByStdName(self,  stdName):
282        '''given a CF standard name return the axis name'''
283        understoodNames={'grid_longitude': 'grid_longitude', 'grid_latitude':'grid_latitude'}
284        if stdName in understoodNames.keys():
285            axisname= understoodNames[stdName]
286        else:
287            axisname =None
288        return axisname
289       
290def main():
291    import cdms # just for testing
292    #sys.stdout = open("/home/dom/svn/trunk/csml/testfiles/gridseries/crsout.txt", "w")
293    cat=CRSCatalogue()
294 
295    #test every crs in this netcdf file
296    crsdictionary={}
297    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc')
298    for var in f.variables:
299        v=f.variables[var]
300        axes=v.getAxisIds()
301        units =[]
302        for axis in axes:
303            ax=f.getAxis(axis)
304            units.append(ax.attributes['units'])
305        crsdictionary[str(axes)]=units
306    n=0
307    #test each definition in dictionary
308       
309    for key in crsdictionary:
310        n=n+1
311        print '\n TEST %d'%n
312        print '*************************************'
313        #make list from string.
314        axs=[]
315        for item in key.split():
316            axs.append(item)
317        axs=eval(key)
318        uns=crsdictionary[key]
319        print 'AXES: %s'%axs
320        print 'UNITS: %s'%uns
321        print '----------------------------'
322        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
323        #print 'CRS: %s'%crs.srsName
324        print 'ORDER:  %s'%axisorder
325        print 'LABELS: %s'%crs.axisLabels
326        print 'CRS AXES: %s'%crs.axes
327        labels=''
328       
329       
330        for axis in axes:
331            labels=labels + axis + ' '
332       
333        #print '\n Pseudo CSML:'
334        #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
335        #print '<axisLabels>%s</axisLabels>'%labels
336        #for axis in axes:
337            #print '<GridOrdinate>'
338            #unit=units[axes.index(axis)]
339            #print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
340            #print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
341            #print '</GridOrdinate>'
342        #print '</GridSeriesDomain>'
343       
344        #print '**************************************'
345   #test getting CRS by name:
346    #print 'getting ndg:crs:xypt by name'
347    #crs=cat.getCRS('ndg:crs:xypt')
348    #print 'returns %s'%crs.srsName
349       
350   
351if __name__=="__main__":
352    main()
Note: See TracBrowser for help on using the repository browser.