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

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

Added new crs (xyh) and integrated into ProfileSeries? code

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.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       
22        # define lon lat pressure time CRS:       
23        crs=CRSystem(srsName='ndg:crs:xypt', axes =['longitude', 'latitude','pressure','time'])
24        crs.lonAxis=0
25        crs.latAxis=1
26        crs.timeAxis=3
27        self.systems['ndg:crs:xypt']=crs
28       
29        # define lon lat height time CRS:       
30        crs=CRSystem(srsName='ndg:crs:xyht', axes =['longitude', 'latitude','height','time'])
31        crs.lonAxis=0
32        crs.latAxis=1
33        crs.timeAxis=3
34        self.systems['ndg:crs:xyht']=crs
35       
36        # define lon lat time CRS:
37        crs=CRSystem(srsName='ndg:crs:xyt', axes =['longitude', 'latitude','time'])
38        self.systems['ndg:crs:xyt']=crs
39        crs.lonAxis=0
40        crs.latAxis=1
41        crs.timeAxis=2
42       
43        # define lon lat height CRS:
44        crs=CRSystem(srsName='ndg:crs:xyh', axes =['longitude', 'latitude','height'])
45        self.systems['ndg:crs:xyh']=crs
46        crs.lonAxis=0
47        crs.latAxis=1
48        crs.timeAxis=None
49       
50       
51        # define height time CRS:
52        crs=CRSystem(srsName='ndg:crs:ht', axes =['height', 'time'])
53        self.systems['ndg:crs:ht']=crs
54        crs.lonAxis=None
55        crs.latAxis=None
56        crs.timeAxis=1
57       
58        #define unknown 1D CRS:
59        crs=CRSystem(srsName='ndg:crs:unknown', axes=['unknown'])
60        self.systems['ndg:crs:unknown']=crs
61       
62        #define unknown 3D CRS:
63        crs=CRSystem(srsName='ndg:crs:unknown2d', axes=['unknown','unknown'])
64        self.systems['ndg:crs:unknown2d']=crs
65       
66        #define unknown 3D CRS:
67        crs=CRSystem(srsName='ndg:crs:unknown3d', axes=['unknown','unknown','unknown'])
68        self.systems['ndg:crs:unknown3d']=crs
69       
70        #define unknown 4D CRS:
71        crs=CRSystem(srsName='ndg:crs:unknown4d', axes=['unknown','unknown','unknown','unknown'])
72        self.systems['ndg:crs:unknown4d']=crs
73       
74           
75    def getCRS(self, srsName):
76        #given the name of a CRS e.g. 'ndg:crs:xypt' return the CRSystem object
77        try:
78            return self.systems[srsName]
79        except:
80            return self.systems['ndg:crs:unknown']
81   
82    def getUnitType(self, unit):
83        unittype='unknown'
84        if string.lower(unit) in ['second', 'seconds', 's', 'mins','minute','minutes','hour','hours','h','hr','hrs','day','days']:
85            unittype='time'
86        elif string.lower(unit)[:10] in ['days since', 'seconds si', 'minutes si', 'hours sinc','months sin', 'years sinc']:
87            unittype='time'
88        elif string.lower(unit) in ['mbar', 'pa','level']:
89            unittype='pressure'
90        elif string.lower(unit) in ['m', 'km']:
91            unittype='height'
92        elif string.lower(unit) in ['degrees_north', 'degrees_south']:
93            unittype='latitude'
94        elif string.lower(unit) in ['degrees_east','degrees_west']:
95            unittype='longitude'
96        return unittype
97   
98    def determineCRS(self, axes=None, units=None, knownCRSAxes=None):
99        '''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.
100        e.g passing in: (axes=['t', 'ht', 'latitude', 'longitude'],units=['days since 1991-09-01 00:00:00', 'm', 'degrees_north', 'degrees_east'])
101        returns self.systems['ndg:crs:xypt'], [3,2,1,0]
102        Assumes the units are ordered to correspond with the axes
103        Alternatively, supply with a list of knownCRSAxes names (e.g. 'latitude', 'longitude', 'height') to get the right crs object       
104        '''
105                       
106        axisorder=[]
107        if knownCRSAxes is not None:
108            axes=knownCRSAxes
109                   
110        if len(axes)==2:
111            crs=self.systems['ndg:crs:unknown2d']
112        if len(axes)==3:
113            crs=self.systems['ndg:crs:unknown3d']
114        elif len(axes)==4:
115            crs=self.systems['ndg:crs:unknown4d']
116        else :
117            crs=self.systems['ndg:crs:unknown']
118       
119        #this can be extended to accomodate more units       
120        if knownCRSAxes is not None:
121            crsMap=knownCRSAxes         
122        else:
123            crsMap=[]
124            for axis in axes:
125                unit=units[axes.index(axis)]
126                unittype=self.getUnitType(unit)           
127                crsMap.append(unittype)
128       
129       
130        #now try and match up crsMap with known crsystems.       
131        #axisorder stores the order of the supplied axes against the known crs.
132        #so if the known crs is 'Lat' 'Lon' 'Time' and the supplied axes are 'Time' 'Lat' 'Lon' then
133        #axisorder would be [2, 0, 1]         
134        for system in self.systems:
135            if len(self.systems[system].axes) != len(crsMap):
136                #wrong dimensionality, skip this crs
137                continue
138            match=0
139            axisorder={}
140            for ax in self.systems[system].axes:
141                if ax in crsMap:
142                    axisorder[ax]=crsMap.index(ax)
143                else:
144                    match = 1
145                    break
146            if match == 0: 
147                crs=self.systems[system]
148                break
149                       
150        return crs, axisorder
151       
152def main():
153    import sys
154    sys.stdout = open("/home/dom/Desktop/SVN/trunk/csml/testfiles/gridseries/crsout.txt", "w")
155    cat=CRSCatalogue()
156 
157    #test every crs in this netcdf file
158    crsdictionary={}
159    f=cdms.open('/home/dom/Desktop/SVN/trunk/csml/testfiles/gridseries/xaaqda@pxs19c1.nc')
160    for var in f.variables:
161        v=f.variables[var]
162        axes=v.getAxisIds()
163        units =[]
164        for axis in axes:
165            ax=f.getAxis(axis)
166            units.append(ax.attributes['units'])
167        crsdictionary[str(axes)]=units
168    n=2
169    #test each definition in dictionary
170    for key in crsdictionary:
171        n=n+1
172        print '\n TEST %d'%n
173        print '*************************************'
174       #make list from string.
175        axs=[]
176        for item in key.split():
177            axs.append(item)
178        uns=crsdictionary[key]
179        print 'AXES: %s'%axs
180        print 'UNITS: %s'%uns
181        crs, axisorder= cat.determineCRS(axes=axs, units=uns)
182        print 'CRS: %s'%crs.srsName
183        print 'ORDER:  %s'%axisorder
184        #print 'axis order =%s'%axisorder
185        labels=''
186       
187        #print axes
188        # ['t', 'ht', 'latitude_1', 'longitude_1']
189        #print crs.axes
190         # ['Lon', 'Lat', 'Pressure', 'Time']
191        #print axisorder
192        # [3, 2, 1, 0]
193       
194        for axis in axes:
195            labels=labels + axis + ' '
196       
197        #print '\n Pseudo CSML:'
198        #print '<GridSeriesDomain srsName=%s srsDimension=%s axisLabels=%s>'%(crs.srsName, crs.srsDimension,crs.axisLabels)
199        #print '<axisLabels>%s</axisLabels>'%labels
200        #for axis in axes:
201            #print '<GridOrdinate>'
202 
203            #unit=units[axes.index(axis)]
204           # print '     <coordAxisLabel>%s</coordAxisLabel>'%crs.axes[axisorder[axes.index(axis)]]
205          #  print '     <gridAxesSpanned>%s</gridAxesSpanned>'%axis
206         #   print '</GridOrdinate>'
207        #print '</GridSeriesDomain>'
208       
209        print '**************************************'
210   #test getting CRS by name:
211    #print 'getting ndg:crs:xypt by name'
212    #crs=cat.getCRS('ndg:crs:xypt')
213    #print 'returns %s'%crs.srsName
214       
215   
216if __name__=="__main__":
217    main()
Note: See TracBrowser for help on using the repository browser.