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

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

fixed bounding box aggregation code

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, knownCRSAxes=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        crsMap=None
238       
239        #stdNames is a list, but may be a list of Nones, in which case set = None.
240        if type(stdNames) is list:
241            namesok=False
242            for item in stdNames:
243                if item is not None:
244                    namesok=True
245            if namesok is False:
246                stdNames=None
247       
248                       
249               
250        if knownCRSAxes is not None:
251            for item in knownCRSAxes:
252                if item is not None:
253                    crsMap=knownCRSAxes
254                    break         
255        if not crsMap:
256            #build a map of axes and units - use std names to help!
257            crsMap=[]
258            for axis in axes:
259                unittype = None
260                if stdNames is not None:
261                    if stdNames[axes.index(axis)]:
262                        unittype=stdNames[axes.index(axis)]
263                else:
264                    unit=units[axes.index(axis)]
265                    unittype=self.getUnitType(unit)         
266                crsMap.append(unittype)             
267        match=0
268        #now try and match up crsMap with known crsystems.       
269        #axisorder stores the order of the supplied axes against the known crs.
270        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
271        #axisorder would be [2, 0, 1]         
272        for system in self.systems:
273            if len(self.systems[system].axes) != len(crsMap):
274                #wrong dimensionality, skip this crs
275                continue
276            match=1
277            axisorder={}
278            for ax in self.systems[system].axes:
279                if ax in crsMap:
280                    axisorder[ax]=crsMap.index(ax)
281                else:
282                    match = 0
283                    break
284            if match == 1: 
285                crs=self.systems[system]             
286                break
287       
288        if match==0:  #if couldn't find a match - generate a pseudo crs     
289            crs, axisorder=self.__generatePseudoCRS(crsMap, axes, units)
290       
291        if crs.units is None:
292            #if no units are defined, set the units to be those of the original data
293            if units is not None:
294                orderedunits=[]
295                for a in crs.axes:
296                    orderedunits.append(units[axisorder[a]])
297                    crs.units=csml.csmllibs.csmlextra.stringify(orderedunits)
298        return crs, axisorder
299   
300    def _identifyAxisByStdName(self,  stdName):
301        '''given a CF standard name return the axis name'''
302        understoodNames={'grid_longitude': 'grid_longitude', 'grid_latitude':'grid_latitude'}
303        if stdName in understoodNames.keys():
304            axisname= understoodNames[stdName]
305        else:
306            axisname =None
307        return axisname
308       
309def main():
310    import cdms # just for testing
311    #sys.stdout = open("/home/dom/svn/trunk/csml/testfiles/gridseries/crsout.txt", "w")
312    cat=CRSCatalogue()
313 
314    #test every crs in this netcdf file
315    crsdictionary={}
316    f=cdms.open('/home/dom/svn/trunk/csml/testfiles/gridseries/xaaqda@pxs09c1.nc')
317    for var in f.variables:
318        v=f.variables[var]
319        axes=v.getAxisIds()
320        units =[]
321        for axis in axes:
322            ax=f.getAxis(axis)
323            units.append(ax.attributes['units'])
324        crsdictionary[str(axes)]=units
325    n=0
326    #test each definition in dictionary
327       
328    for key in crsdictionary:
329        n=n+1
330        print '\n TEST %d'%n
331        print '*************************************'
332        #make list from string.
333        axs=[]
334        for item in key.split():
335            axs.append(item)
336        axs=eval(key)
337        uns=crsdictionary[key]
338        print 'AXES: %s'%axs
339        print 'UNITS: %s'%uns
340        print '----------------------------'
341        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
342        #print 'CRS: %s'%crs.srsName
343        print 'ORDER:  %s'%axisorder
344        print 'LABELS: %s'%crs.axisLabels
345        print 'CRS AXES: %s'%crs.axes
346        labels=''
347       
348       
349        for axis in axes:
350            labels=labels + axis + ' '
351       
352        #print '\n Pseudo CSML:'
353        #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
354        #print '<axisLabels>%s</axisLabels>'%labels
355        #for axis in axes:
356            #print '<GridOrdinate>'
357            #unit=units[axes.index(axis)]
358            #print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
359            #print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
360            #print '</GridOrdinate>'
361        #print '</GridSeriesDomain>'
362       
363        #print '**************************************'
364   #test getting CRS by name:
365    #print 'getting ndg:crs:xypt by name'
366    #crs=cat.getCRS('ndg:crs:xypt')
367    #print 'returns %s'%crs.srsName
368       
369   
370if __name__=="__main__":
371    main()
Note: See TracBrowser for help on using the repository browser.