Module sia_model
[hide private]
[frames] | no frames]

Source Code for Module sia_model

  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  # standard libraries 
  7  import sys 
  8  import os 
  9  import copy 
 10   
 11  # numerics 
 12  import numpy as N 
 13  import pylab as P 
 14  ## this is the old numerics module, outdate, but we need it because of 
 15  ## the netCDF libaray we use 
 16  import Numeric 
 17   
 18  # our modules 
 19  import helper_module 
 20  sys.path.append('fortran') 
 21  import kees_simple_ice_model 
 22   
23 -class Structure(object):
24 """ 25 An empty class serving as a container. Similar to a derived type 26 in Fortran. 27 """ 28 pass
29
30 -class SIA_model(object):
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 #: dict holding the L{Dimension} instances 45 self.dimensions = {} 46 #: dict holding the L{Variable} instances 47 self.variables = {} 48 #: dict holding the global attributes 49 self.global_attrs = {} 50 51 if conf_or_netcdf_file[-3:]=='.nc': 52 # if a .nc file read the netCDF file 53 self.read_netCDF() 54 elif conf_or_netcdf_file[-4:]=='.CDL': 55 # if a .CDL file read the CDL file 56 self.read_model_config() 57 self.read_model_input()
58 59 60
61 - def read_model_config(self):
62 """ 63 Read the configuration file of the model. 64 """ 65 tmp_file = 'tmp/config.nc' 66 # first remove old files from the tmp folder 67 os.system('rm '+ tmp_file) 68 # turn the CDL file into a .nc file 69 helper_module.CDL2nc(self.input_file, tmp_file) 70 # read it 71 config_file = helper_module.open_netCDF_into_python(tmp_file) 72 73 ## now comes the hard bit: fill the class with 74 ## the specified data 75 76 # global attrs: just set them to what they are in the config 77 # file 78 self.global_attrs = self.fix_Scientific_bug(config_file.__dict__) 79 80 # read the input files and put their data into the appropriate 81 # place 82 for key in self.global_attrs: 83 ## make input variables 84 if key.startswith('input_file'): 85 input_file = helper_module.open_netCDF_into_python( 86 self.global_attrs[key]) 87 # here we assume that the dimension (X) is the same for all 88 # input and of the same length. Thus set the dimension 89 # to that length: 90 if self.dimensions=={}: 91 self.createDimension( 92 'X', input_file.dimensions['X']) 93 # and the corresponding coordinate variable 94 X = input_file.variables['X'] 95 # make into numpy type, see note below: 96 dict_ = self.fix_Scientific_bug(X.__dict__) 97 self.createVariable('X', X.typecode(), 98 dimensions=('X',), attrs=dict_) 99 # and set the data 100 self.variables['X'].data = N.array(X[:]) 101 # note: X[:] is not a numpy array, but a older 102 # numeric array, thus I turn it into a numpy 103 # array. 104 105 # add the variable 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 # and set the data 114 self.variables[kk].data = N.array(var[:]) 115 116 # note that we just discarded the global attributs of 117 # the input files 118 # finally, close the config file 119 config_file.close() 120 121 ## make the time dimension and variable: 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 ## make the output variable 130 thick_attrs = {'units': 'm', 'standard_name': 'land_ice_thickness'} 131 self.createVariable('thick', 'd', ('X', 'time'), thick_attrs) 132 133 # for the fortran module to know when ot output we make t_out_ind: 134 # the indices of the time steps (dt) when to output 135 self.time_steps_between_output = N.round(self.global_attrs['t_output_dt']/ 136 self.global_attrs['dt']) 137 #: number of integeration steps 138 self.t_steps = N.round(self.global_attrs['t_final']/ 139 self.global_attrs['dt'])
140
141 - def read_model_input(self):
142 """Reads all the input files: 143 - bed 144 - initial thickness 145 - mass balance distribution 146 """ 147 pass
148
149 - def read_netCDF(self):
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
157 - def write_netCDF(self):
158 """ 159 Write all the calcualted data and also metadata to a netCDF 160 file. 161 """ 162 # open the file specified in the config file 163 netCDFfile = helper_module.open_netCDF_into_python( 164 self.global_attrs['output_filename'], 'w') 165 # now follow U{http://gfesuite.noaa.gov/developer/netCDFPythonInterface.html} 166 167 ## global attributs 168 for kk in self.global_attrs: 169 setattr(netCDFfile, kk, 5) 170 171 ## first the dimensions 172 for kk in self.dimensions: 173 netCDFfile.createDimension(kk, self.dimensions[kk].size) 174 ## now the variables 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 # fill in the data 181 print kk 182 print netCDFvar.dimensions 183 print var.name, var.data.shape 184 netCDFvar[:] = var.data 185 ## write to disk 186 netCDFfile.sync()
187
188 - def run_model(self):
189 """ 190 Runs the model according to the setup gathered from 191 L{read_model_config} and L{read_input}. 192 """ 193 ##thick_out = simple_ice_model(dt,t_steps,n_t_num,dx,rho,g_grav,a_rate,n_glen,xx,mass_b,bed,thick,thickr,thickl,[grid]) 194 195 ## in fortran 196 # subroutine simple_ice_model( dt, t_steps, n_t_num, dx, & ! numerics: 197 # ! time step, total number of time steps, time steps between outputs, spatial discretisation step 198 # rho, g_grav, A_rate, n_glen, & ! physical constants 199 # xx, mass_b, bed, & ! input on grid: the grid, mass balance, bed elevation 200 # thick, & !IC on grid: initial ice thickness 201 # thickr, thickl, & ! BC: thickness at left and right 202 # thick_out) ! output: the thickness at the specified time steps 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
218 - def plot_thickness(self, time):
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 # get what index is the solution at 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
247 - def plot_IC(self):
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 # helper methods xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
276 - def createDimension(self, dimName, size):
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
308 - def fix_Scientific_bug(self, input_dict):
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
326 -class Variable(object):
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 #: the variable name 355 self.name = varname 356 #: the datatype 357 self.datatype = datatype 358 #: netCDF CF-1.4 attributes 359 self.attrs = attrs 360 361 self.dimensions = dimensions 362 #: the acutal data in a numpy array (this could probably be done better) 363 self.data = N.array([])
364
365 - def get_axis_label(self):
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
374 -class Dimension(object):
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
385 - def __init__(self, dimname, size):
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