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

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

improved determineCRS method

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