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