1 | '''#getSubset.py contains generic subsetting code that can be reused for multiple features |
---|
2 | Dominic Lowe, CCLRC, 31/01/2007''' |
---|
3 | |
---|
4 | import csml |
---|
5 | import csmlutils |
---|
6 | |
---|
7 | |
---|
8 | def subsetDomain(timeaxis,times, domain,**kwargs): |
---|
9 | '''takes the domain and returns a subset according to the keyword selection criteria |
---|
10 | time = name of time dimension |
---|
11 | times = string of times |
---|
12 | ''' |
---|
13 | strTimes=times |
---|
14 | subsetDomain={} |
---|
15 | totalArraySize=1 |
---|
16 | for key in domain.keys(): |
---|
17 | straxisValues='' |
---|
18 | if key in kwargs: |
---|
19 | if key ==timeaxis: |
---|
20 | straxisValues=strTimes |
---|
21 | arraySize=len(strTimes.split()) |
---|
22 | elif kwargs[key][0] < kwargs[key][1]: |
---|
23 | for val in domain[key]: |
---|
24 | if val is not '': |
---|
25 | if float(val) >= kwargs[key][0]: |
---|
26 | if float(val) <= kwargs[key] [1]: |
---|
27 | straxisValues=straxisValues+ str(val) + ',' |
---|
28 | straxisValues=straxisValues[:-1] |
---|
29 | else:#this if deals with selections such as longitude (330,30) where the lower limit is 'greater' than the upper limit in a mathematical sense. |
---|
30 | for val in domain[key]: |
---|
31 | if val is not '': |
---|
32 | if val >= kwargs[key][0]: |
---|
33 | straxisValues=straxisValues+ str(val) + ',' |
---|
34 | for val in domain[key]: |
---|
35 | if val is not '': |
---|
36 | if val <= kwargs[key] [1]: |
---|
37 | straxisValues=straxisValues+ str(val) + ',' |
---|
38 | straxisValues=straxisValues[:-1] |
---|
39 | else: # this dimension has not been subsetted at all - CHECK THIS |
---|
40 | for val in domain[key]: |
---|
41 | if val is not '': |
---|
42 | straxisValues=straxisValues+ str(val) + ',' |
---|
43 | straxisValues=straxisValues[:-1] |
---|
44 | |
---|
45 | subsetDomain[key]=straxisValues |
---|
46 | if key != timeaxis: |
---|
47 | arraySize=len(subsetDomain[key].split(',')) |
---|
48 | else: |
---|
49 | arraySize=len(subsetDomain[key].split()) |
---|
50 | totalArraySize = totalArraySize * arraySize |
---|
51 | return subsetDomain, totalArraySize |
---|
52 | |
---|
53 | def _setCalendar(feature, timeName): |
---|
54 | '''sets the cdtime calendar needed for conversion from CSML to NetCDF''' |
---|
55 | calset=False |
---|
56 | for gridOrd in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates: |
---|
57 | if gridOrd.coordAxisLabel.CONTENT==timeName: |
---|
58 | try: |
---|
59 | caltype=gridOrd.coordAxisValues.frame.split(':',1)[0] |
---|
60 | calunits=gridOrd.coordAxisValues.frame.split(':',1)[1] |
---|
61 | csml.csmllibs.csmltime.setcdtimeCalendar(caltype) |
---|
62 | calset=True |
---|
63 | except:pass |
---|
64 | if calset!=True: |
---|
65 | csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar) |
---|
66 | try: |
---|
67 | caltype=gridOrd.coordAxisValues.frame.split(':',1)[0] |
---|
68 | calunits=gridOrd.coordAxisValues.frame.split(':',1)[1] |
---|
69 | csml.csmllibs.csmltime.setcdtimeCalendar(caltype) |
---|
70 | except: |
---|
71 | csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar) |
---|
72 | return calunits, caltype |
---|
73 | |
---|
74 | |
---|
75 | def _getTimes(timeSelection, timeName, domain): |
---|
76 | |
---|
77 | #currently supporting domain subsetting only by CRS name |
---|
78 | #(but should be easy to extend later) |
---|
79 | times=timeSelection |
---|
80 | |
---|
81 | tone=csml.csmllibs.csmltime.getCDtime(times[0]) |
---|
82 | if len(times) == 2: |
---|
83 | #then this is a range of times (t1, tn), and not a list |
---|
84 | try: |
---|
85 | tone=csml.csmllibs.csmltime.getCDtime(times[0]) |
---|
86 | except: |
---|
87 | tone=times[0] |
---|
88 | try: |
---|
89 | ttwo=csml.csmllibs.csmltime.getCDtime(times[1]) |
---|
90 | except: |
---|
91 | ttwo=times[1] |
---|
92 | times=[] |
---|
93 | for time in domain[timeName]: |
---|
94 | timeok=csml.csmllibs.csmltime.compareTimes(tone,time,ttwo) |
---|
95 | if timeok ==1: |
---|
96 | times.append(time) |
---|
97 | return times |
---|
98 | |
---|
99 | def _getTimeToFileRatio(feature,domain, timeName): |
---|
100 | if hasattr(feature.value.rangeSet, 'valueArray'): |
---|
101 | if hasattr(feature.value.rangeSet.valueArray.valueComponent.quantityList, '__insertedExtract'): |
---|
102 | numFiles= len( csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components)[0].fileList.fileNames.CONTENT.split()) |
---|
103 | return len(domain[timeName])/numFiles |
---|
104 | |
---|
105 | def getCoordTransformTable(domainSubset, crs): |
---|
106 | cTT=csml.parser.GridCoordinatesTable() |
---|
107 | ords =[] |
---|
108 | for key in domainSubset.keys(): |
---|
109 | go=csml.parser.GridOrdinateDescription() |
---|
110 | go.coordAxisLabel=csml.parser.csString(key) |
---|
111 | go.gridAxesSpanned=csml.parser.csString(key) |
---|
112 | go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList() |
---|
113 | if key==crs.axes[crs.timeAxis]: |
---|
114 | go.coordAxisValues.timePositionList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder |
---|
115 | else: |
---|
116 | go.coordAxisValues.coordinateList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder |
---|
117 | ords.append(go) #go needs a few more properties setting |
---|
118 | cTT.gridOrdinates=ords |
---|
119 | return cTT |
---|
120 | |
---|
121 | |
---|
122 | |
---|
123 | def getTheData(feature, selection, times,timeName): |
---|
124 | domain = feature.domain |
---|
125 | value=feature.value |
---|
126 | files=[] |
---|
127 | strTimes='' |
---|
128 | fulldata=[] |
---|
129 | #get the ratio of times to files |
---|
130 | timeToFileRatio = _getTimeToFileRatio(feature, domain, timeName) |
---|
131 | |
---|
132 | |
---|
133 | #list to keep track of files that have already been fetched. eg. if multiple times are in a single file only need to get data from that file once... |
---|
134 | filesFetched=[] |
---|
135 | |
---|
136 | |
---|
137 | for time in times: |
---|
138 | listPosition=domain[timeName].index(time) |
---|
139 | strTimes= strTimes + ' ' + time |
---|
140 | for comp in csmlutils.listify(value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): |
---|
141 | filePos=int(float(listPosition)/timeToFileRatio) |
---|
142 | if filePos in filesFetched: |
---|
143 | continue #already got data from this file, try next time |
---|
144 | data=comp.getData(fileposition=filePos, **selection) |
---|
145 | files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name |
---|
146 | if fulldata ==[]: |
---|
147 | fulldata = data.tolist() |
---|
148 | else: |
---|
149 | for item in data.tolist(): |
---|
150 | fulldata.append(item) |
---|
151 | filesFetched.append(filePos) |
---|
152 | axisorder = data.getAxisIds() #will need later! |
---|
153 | |
---|
154 | return strTimes, axisorder, fulldata |
---|
155 | |
---|
156 | def genericSubset(feature, csmlpath, ncpath, domain, kwargs): |
---|
157 | if ncpath is not None: |
---|
158 | pathToSubsetNetCDF=ncpath |
---|
159 | else: |
---|
160 | pathToSubsetNetCDF='temp.nc' |
---|
161 | |
---|
162 | #deal with longitude requests |
---|
163 | #if the request is in -ve,+ve eg (-30,30) but the data is in (0,360) need to handle this by changing the args. |
---|
164 | kwargs=csmlutils.fixLongitude(domain,kwargs) |
---|
165 | |
---|
166 | |
---|
167 | #get the name of the time axis in the coordinate reference system |
---|
168 | cat=csml.csmllibs.csmlcrs.CRSCatalogue() |
---|
169 | crs=cat.getCRS(feature.value.gridSeriesDomain.srsName) |
---|
170 | timeName=crs.axes[crs.timeAxis] |
---|
171 | #set the reference system for the time axis |
---|
172 | calunits, caltype=_setCalendar(feature, timeName) |
---|
173 | |
---|
174 | #selection by time - get explcit time values, not just a range mint, maxt |
---|
175 | timeSelection=kwargs[timeName] |
---|
176 | times=_getTimes(timeSelection, timeName,domain) |
---|
177 | |
---|
178 | |
---|
179 | return pathToSubsetNetCDF, kwargs, timeName, calunits, caltype, times |
---|