1 """
2 This is a sample class for setting up, running and output of a model.
3 In this instance a shallow ice sheet model.
4 """
5
6
7 import sys
8 import os
9 import copy
10
11
12 import numpy as N
13 import pylab as P
14
15
16 import Numeric
17
18
19 import helper_module
20 sys.path.append('fortran')
21 import kees_simple_ice_model
22
24 """
25 An empty class serving as a container. Similar to a derived type
26 in Fortran.
27 """
28 pass
29
31 """
32 The class doing all the model setup, running, output and some
33 plotting.
34 """
35 - def __init__(self, conf_or_netcdf_file):
36 """The standard initialisation function.
37
38 @type conf_or_netcdf_file: string
39 @param conf_or_netcdf_file: A config file or a former netCDF
40 output file.
41 """
42 self.input_file = conf_or_netcdf_file
43
44
45 self.dimensions = {}
46
47 self.variables = {}
48
49 self.global_attrs = {}
50
51 if conf_or_netcdf_file[-3:]=='.nc':
52
53 self.read_netCDF()
54 elif conf_or_netcdf_file[-4:]=='.CDL':
55
56 self.read_model_config()
57 self.read_model_input()
58
59
60
62 """
63 Read the configuration file of the model.
64 """
65 tmp_file = 'tmp/config.nc'
66
67 os.system('rm '+ tmp_file)
68
69 helper_module.CDL2nc(self.input_file, tmp_file)
70
71 config_file = helper_module.open_netCDF_into_python(tmp_file)
72
73
74
75
76
77
78 self.global_attrs = self.fix_Scientific_bug(config_file.__dict__)
79
80
81
82 for key in self.global_attrs:
83
84 if key.startswith('input_file'):
85 input_file = helper_module.open_netCDF_into_python(
86 self.global_attrs[key])
87
88
89
90 if self.dimensions=={}:
91 self.createDimension(
92 'X', input_file.dimensions['X'])
93
94 X = input_file.variables['X']
95
96 dict_ = self.fix_Scientific_bug(X.__dict__)
97 self.createVariable('X', X.typecode(),
98 dimensions=('X',), attrs=dict_)
99
100 self.variables['X'].data = N.array(X[:])
101
102
103
104
105
106 for kk in input_file.variables.keys():
107 if not kk=='X':
108 var = input_file.variables[kk]
109 dict_ = self.fix_Scientific_bug(var.__dict__)
110 self.createVariable(kk, var.typecode(),
111 dimensions=('X',),
112 attrs=dict_)
113
114 self.variables[kk].data = N.array(var[:])
115
116
117
118
119 config_file.close()
120
121
122 time = N.arange(self.global_attrs['t_start'], self.global_attrs['t_final'],
123 self.global_attrs['t_output_dt'])
124 self.createDimension('time', len(time))
125 time_attrs = {'units': self.global_attrs['t_units'], 'standard_name': 'time'}
126 time_var = self.createVariable('time', 'd', ('time',), time_attrs)
127 time_var.data = time
128
129
130 thick_attrs = {'units': 'm', 'standard_name': 'land_ice_thickness'}
131 self.createVariable('thick', 'd', ('X', 'time'), thick_attrs)
132
133
134
135 self.time_steps_between_output = N.round(self.global_attrs['t_output_dt']/
136 self.global_attrs['dt'])
137
138 self.t_steps = N.round(self.global_attrs['t_final']/
139 self.global_attrs['dt'])
140
148
150 """
151 Reads a netCDF file (as saved by a previous model run) and
152 puts everything into place, such that the model run can be
153 contiued.
154 """
155 pass
156
158 """
159 Write all the calcualted data and also metadata to a netCDF
160 file.
161 """
162
163 netCDFfile = helper_module.open_netCDF_into_python(
164 self.global_attrs['output_filename'], 'w')
165
166
167
168 for kk in self.global_attrs:
169 setattr(netCDFfile, kk, 5)
170
171
172 for kk in self.dimensions:
173 netCDFfile.createDimension(kk, self.dimensions[kk].size)
174
175 for kk in self.variables:
176 var = self.variables[kk]
177 for kkk in var.attrs:
178 setattr(netCDFfile, kkk, 5)
179 netCDFvar = netCDFfile.createVariable(kk, var.datatype, var.dimensions)
180
181 print kk
182 print netCDFvar.dimensions
183 print var.name, var.data.shape
184 netCDFvar[:] = var.data
185
186 netCDFfile.sync()
187
189 """
190 Runs the model according to the setup gathered from
191 L{read_model_config} and L{read_input}.
192 """
193
194
195
196
197
198
199
200
201
202
203
204 X = self.variables['X']
205 land_ice_copy = copy.deepcopy(self.variables['land_ice_thickness'].data)
206 dx = X.data[1]-X.data[0]
207 time = self.variables['time']
208 print self.time_steps_between_output
209 output_thickness = kees_simple_ice_model.simple_ice_model(
210 self.global_attrs['dt'], self.t_steps, self.time_steps_between_output, dx,
211 self.global_attrs['rho_i'], self.global_attrs['g_grav'], self.global_attrs['A_rate'], self.global_attrs['n_glen'],
212 X.data, self.variables['land_ice_surface_specific_mass_balance'].data, self.variables['bedrock_altitude'].data,
213 land_ice_copy,
214 self.global_attrs['BC_thickl'], self.global_attrs['BC_thickr'], len(X.data))
215
216 self.variables['thick'].data = output_thickness
217
219 """
220 Plot the ice thickness at a particular time on top of the bed.
221
222 @type time: float
223 @param time: The time when to plot the solution. If <0 then
224 the number is used as index (from the last).
225 """
226
227 if time<0:
228 ind = time
229 time = self.variables['time'].data[-1]
230 else:
231 ind = N.where(time>self.variables['time'].data)[0][-1]
232
233 fig = P.figure()
234 P.plot(self.variables['X'].data,
235 self.variables['bedrock_altitude'].data +
236 self.variables['thick'].data[:,ind].T, 'b')
237 P.hold(True)
238 P.plot(self.variables['X'].data,
239 self.variables['bedrock_altitude'].data, 'r',
240 label=str(self.variables['bedrock_altitude'].attrs['standard_name']))
241 P.xlabel(self.variables['time'].get_axis_label())
242 P.ylabel('Elevation (%s)' % self.variables['bedrock_altitude'].attrs['units'])
243 P.title('Solution at time %.0f (%s)'
244 %(time, self.variables['time'].attrs['units']))
245
246
248 """
249 Plot the ice thickness at t0, together with the bed and mass
250 balance distribution.
251 """
252 fig = P.figure()
253 ax1 = P.subplot(2,1,1)
254 P.plot(self.variables['X'].data,
255 self.variables['land_ice_thickness'].data +
256 self.variables['bedrock_altitude'].data , 'b',
257 label=str(self.variables['land_ice_thickness'].attrs['standard_name']))
258 P.hold(True)
259 P.plot(self.variables['X'].data,
260 self.variables['bedrock_altitude'].data, 'r',
261 label=str(self.variables['bedrock_altitude'].attrs['standard_name']))
262 P.xlabel(self.variables['X'].get_axis_label())
263 P.ylabel('Elevation (%s)' % self.variables['bedrock_altitude'].attrs['units'])
264 P.legend()
265 P.title('IC and other input')
266
267 ax1 = P.subplot(2,1,2)
268 P.plot(self.variables['X'].data,
269 self.variables['land_ice_surface_specific_mass_balance'].data)
270 P.xlabel(self.variables['X'].get_axis_label())
271 P.ylabel(self.variables['land_ice_surface_specific_mass_balance'].get_axis_label())
272
273
274
275
277 """
278 Add a L{Dimension} instance.
279
280 @type dimName: string
281 @param dimName: The name of the dimension
282
283 @type size: integer
284 @param size: the length of the dimension
285 """
286 self.dimensions[dimName] = Dimension(dimName, size)
287 return self.dimensions[dimName]
288
289 - def createVariable(self, varname, datatype='d', dimensions=(), attrs={}):
290 """
291 Add a L{Dimension} instance.
292
293 @type varname: string
294 @param varname: The name of the variable
295
296 @type dimensions: a touple of strings
297 @param dimensions: what dimensions are associated with this
298 variable, if () then its scalar.
299
300 @type attrs: dict
301 @param attrs: a dictionary holding all the attirbutes of the
302 variable.
303 """
304 self.variables[varname] = Variable(varname, datatype
305 , dimensions, attrs)
306 return self.variables[varname]
307
309 """
310 The C{Scientific.IO.NetCDF.NetCDFFile} retruns not numpy
311 arrays but other one of a now outdated package numeric. Thus
312 all needs to be turned into numpy arrays.
313
314 @type input_dict: dict
315 @param input_dict: A dictionary containing numeric arrays
316
317 @return: A dictionary containing numpy arrays
318 """
319 dict_ = {}
320 for kk in input_dict:
321 dict_[kk] = N.array(input_dict[kk])
322 return dict_
323
324
325
327 """
328 A class to hold variable data and attributes, similar to the
329 C{Scientific.io.netcdf} variables.
330 """
331
332 - def __init__(self,varname, datatype='d', dimensions=(), attrs={}):
333 """
334 @type varname: string
335 @param varname: Name of the variable.
336
337 @type datatype: string
338 @param datatype: One of:
339 - 'f4' (32-bit floating point)
340 - 'd' (64-bit floating point)
341 - 'i4' (32-bit signed integer)
342 - 'i2' (16-bit signed integer)
343 - 'i1' (8-bit signed integer)
344 - 'S1' (single-character string)
345
346 @type dimensions: a touple of strings
347 @param dimensions: what dimensions are associated with this
348 variable, if () then its scalar.
349
350 @type attrs: dict
351 @param attrs: a dictionary holding all the attirbutes of the
352 variable.
353 """
354
355 self.name = varname
356
357 self.datatype = datatype
358
359 self.attrs = attrs
360
361 self.dimensions = dimensions
362
363 self.data = N.array([])
364
366 """
367 Returns a string to label a axis in a plot.
368 """
369 lab = '%s (%s)' %(self.attrs['standard_name'],
370 self.attrs['units'])
371 return lab
372
373
375 """
376 A class to hold a dimension.
377
378 Note that each (netCDF) variable has a number of dimensions
379 associated with it, i.e. a touple containing instances of this
380 class. The dimension should have a closely associated variable,
381 I{coordinate variable} which needs to have the same name as the
382 dimension.
383 """
384
386 """
387 @type dimname: string
388 @param dimname: The standart name according CF 1.4 convention.
389
390 @type size: int
391 @param size: The length of the coordinate variable.
392 """
393 self.name = dimname
394 self.size = size
395