1 | """Create mesh and time boundary for the Okushiri island validation |
---|
2 | """ |
---|
3 | |
---|
4 | |
---|
5 | from Numeric import array, zeros, Float, allclose |
---|
6 | |
---|
7 | from anuga.pmesh.mesh import * |
---|
8 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
9 | from anuga.geospatial_data import Geospatial_data |
---|
10 | |
---|
11 | import project_truescale |
---|
12 | |
---|
13 | def prepare_bathymetry(filename): |
---|
14 | """Convert benchmark 2 bathymetry to NetCDF pts file. |
---|
15 | This is a 'throw-away' code taylor made for files like |
---|
16 | 'okushiri_truescale_bathymetry.txt' from the LWRU2 benchmark |
---|
17 | """ |
---|
18 | |
---|
19 | print 'Creating', filename |
---|
20 | |
---|
21 | # Read the ascii (.txt) version of this file, |
---|
22 | # make it comma separated and invert the bathymetry |
---|
23 | # (Below mean sea level should be negative) |
---|
24 | infile = open(filename[:-4] + '.txt') |
---|
25 | |
---|
26 | points = [] |
---|
27 | attribute = [] |
---|
28 | for line in infile.readlines()[1:]: #Skip first line (the header) |
---|
29 | fields = line.strip().split() |
---|
30 | |
---|
31 | x = float(fields[0]) |
---|
32 | y = float(fields[1]) |
---|
33 | z = -float(fields[2]) # Bathymetry is inverted in original file |
---|
34 | |
---|
35 | points.append([x,y]) |
---|
36 | attribute.append(z) |
---|
37 | infile.close() |
---|
38 | |
---|
39 | # Convert to geospatial data and store as NetCDF |
---|
40 | G = Geospatial_data(data_points=points, |
---|
41 | attributes=attribute) |
---|
42 | G.export_points_file(filename) |
---|
43 | |
---|
44 | |
---|
45 | |
---|
46 | def prepare_timeboundary(filename): |
---|
47 | """Convert benchmark 2 time series to NetCDF tms file. |
---|
48 | This is a 'throw-away' code taylor made for files like |
---|
49 | 'okushiri_truescale_input.txt' from the LWRU2 benchmark |
---|
50 | """ |
---|
51 | |
---|
52 | from Scientific.IO.NetCDF import NetCDFFile |
---|
53 | from Numeric import array |
---|
54 | |
---|
55 | |
---|
56 | print 'Creating', filename |
---|
57 | |
---|
58 | # Read the ascii (.txt) version of this file |
---|
59 | fid = open(filename[:-4] + '.txt') |
---|
60 | |
---|
61 | # Skip first line |
---|
62 | line = fid.readline() |
---|
63 | |
---|
64 | # Read remaining lines |
---|
65 | lines = fid.readlines() |
---|
66 | fid.close() |
---|
67 | |
---|
68 | |
---|
69 | N = len(lines) |
---|
70 | T = zeros(N, Float) #Time |
---|
71 | Q = zeros(N, Float) #Values |
---|
72 | |
---|
73 | for i, line in enumerate(lines): |
---|
74 | fields = line.split() |
---|
75 | |
---|
76 | T[i] = float(fields[0]) |
---|
77 | Q[i] = float(fields[1]) |
---|
78 | |
---|
79 | |
---|
80 | # Create tms NetCDF file |
---|
81 | |
---|
82 | fid = NetCDFFile(filename, 'w') |
---|
83 | fid.institution = 'Geoscience Australia' |
---|
84 | fid.description = 'Input wave for truescale okushiri wavetank experiment' |
---|
85 | fid.starttime = 0.0 |
---|
86 | fid.createDimension('number_of_timesteps', len(T)) |
---|
87 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
88 | fid.variables['time'][:] = T |
---|
89 | |
---|
90 | fid.createVariable('stage', Float, ('number_of_timesteps',)) |
---|
91 | fid.variables['stage'][:] = Q[:] |
---|
92 | |
---|
93 | fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) |
---|
94 | fid.variables['xmomentum'][:] = 0.0 |
---|
95 | |
---|
96 | fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) |
---|
97 | fid.variables['ymomentum'][:] = 0.0 |
---|
98 | |
---|
99 | fid.close() |
---|
100 | |
---|
101 | |
---|
102 | #------------------------------------------------------------- |
---|
103 | if __name__ == "__main__": |
---|
104 | |
---|
105 | |
---|
106 | # Prepare bathymetry |
---|
107 | prepare_bathymetry(project_truescale.bathymetry_filename) |
---|
108 | |
---|
109 | # Prepare time boundary |
---|
110 | prepare_timeboundary(project_truescale.boundary_filename) |
---|
111 | |
---|
112 | |
---|
113 | |
---|