source: development/stochastic_study/create_realisations.py @ 3295

Last change on this file since 3295 was 3035, checked in by ole, 19 years ago

Working with Suresh at ACFR

File size: 2.9 KB
Line 
1"""Create several realisations of bathymetries for use with run_model.
2
3This script will create a number of random samples of the input bathymetry and
4then fit all of them to arrays that can be used in set_quantity inside run_model.
5
6The reason for this is that doing all fits in one sweep saves computation time
7"""
8
9from RandomArray import normal, seed
10from Numeric import concatenate, NewAxis
11import cPickle
12from load_mesh.loadASCII import import_points_file, import_mesh_file
13from pyvolution.least_squares import fit_to_mesh
14import project
15
16bathymetry_filename = project.bathymetry_filename
17mesh_filename = project.working_dir + project.mesh_filename
18
19
20# Read mesh
21mesh_dict = import_mesh_file(mesh_filename)
22vertex_coordinates = mesh_dict['vertices']
23triangles = mesh_dict['triangles']
24
25# Read basic bathymetry from NetCDF .pts file
26print 'Reading bathymetry'
27base_bathymetry = import_points_file(bathymetry_filename)
28
29
30# Extract relevant Numeric arrays
31print 'Extract bathymetry'
32point_coordinates = base_bathymetry['pointlist']
33elevation = base_bathymetry['attributelist']['elevation']
34N = len(elevation)
35print 'read', N, 'points:' + 'z in [%f, %f]' %(min(elevation),
36                                               max(elevation))
37
38# Create different bathymetries based on sampling (Suresh to fill this in)
39print 'Create different realisations'
40Z = None
41block_number = 0
42seed(13,17)
43for n in range(project.number_of_realisations):
44
45    print 'Creating sample #%d' %n
46    # nth realisation of bathymetry
47
48    if n == 0:
49        z = elevation # Keep the first realisation as original
50    else:   
51        z = elevation + normal(project.mean,
52                               project.std_dev,
53                               N)
54
55    #Concatenate into array with each realisation a column
56    if Z is None:
57        Z = z[:, NewAxis] # The first realisation of this block
58    else:   
59        Z = concatenate((Z, z[:, NewAxis]), axis=1)
60
61    if (n+1)%project.blocksize == 0 or n == project.number_of_realisations-1:
62        #Fit all of them to computational mesh
63        print 'fit block %d to the mesh' %block_number
64        V = fit_to_mesh(vertex_coordinates,
65                        triangles,
66                        point_coordinates,
67                        Z,
68                        alpha = 0.0001,
69                        verbose = True,
70                        expand_search = True,
71                        precrop = True,
72                        use_cache = False)
73
74        Z = None
75
76        #Store this block for use with model
77        print 'Store fitted quantities in block %d for use with model'\
78              %block_number
79
80        filename = project.working_dir + project.basename +\
81                   '_block%d.pck' %block_number
82
83        block_number += 1   
84       
85       
86        fid = open(filename, 'wb')
87        cPickle.dump(V, fid)
88        fid.close()
89       
Note: See TracBrowser for help on using the repository browser.