[2424] | 1 | """Read in tms file, interpolate at specified locations (ch 5,7,9) and compare |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | import sys |
---|
| 5 | from os import sep |
---|
| 6 | from caching import cache |
---|
[2425] | 7 | import project |
---|
[2424] | 8 | |
---|
| 9 | sys.path.append('..'+sep+'..'+sep) |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | gauges = [[0.000, 1.696]] #Boundary |
---|
| 13 | gauges += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9 |
---|
| 14 | |
---|
| 15 | gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9'] |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | finaltime = 22.5 |
---|
| 19 | timestep = 0.05 |
---|
| 20 | |
---|
| 21 | #Read reference data |
---|
| 22 | |
---|
| 23 | #Input wave |
---|
[2425] | 24 | filename = 'input_wave.tms' |
---|
[2424] | 25 | print 'Reading', filename |
---|
| 26 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 27 | from Numeric import allclose |
---|
| 28 | fid = NetCDFFile(filename, 'r')# |
---|
| 29 | |
---|
| 30 | input_time = fid.variables['time'][:] |
---|
| 31 | input_stage = fid.variables['stage'][:] |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | #gauges |
---|
| 35 | reference_time = [] |
---|
| 36 | ch5 = [] |
---|
| 37 | ch7 = [] |
---|
| 38 | ch9 = [] |
---|
| 39 | filename = 'output_ch5-7-9.txt' |
---|
| 40 | fid = open(filename) |
---|
| 41 | lines = fid.readlines() |
---|
| 42 | fid.close() |
---|
| 43 | for i, line in enumerate(lines[1:]): |
---|
| 44 | if i == len(input_time): break |
---|
| 45 | |
---|
| 46 | fields = line.split() |
---|
| 47 | |
---|
| 48 | reference_time.append(float(fields[0])) |
---|
| 49 | ch5.append(float(fields[1])/100) #cm2m |
---|
| 50 | ch7.append(float(fields[2])/100) #cm2m |
---|
| 51 | ch9.append(float(fields[3])/100) #cm2m |
---|
| 52 | |
---|
| 53 | |
---|
[2610] | 54 | from pyvolution.util import file_function |
---|
| 55 | from utilities.numerical_tools import ensure_numeric |
---|
[2424] | 56 | gauge_values = [ensure_numeric(input_stage), |
---|
| 57 | ensure_numeric(ch5), |
---|
| 58 | ensure_numeric(ch7), |
---|
| 59 | ensure_numeric(ch9)] #Reference values |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | #Read model output |
---|
[2979] | 64 | #filename = project.basename + '_original.sww' |
---|
| 65 | filename = project.basename + '.sww' |
---|
[2424] | 66 | |
---|
[2425] | 67 | f = cache(file_function, filename, |
---|
[2424] | 68 | {'quantities': 'stage', |
---|
| 69 | 'interpolation_points': gauges, |
---|
| 70 | 'verbose': True}, |
---|
| 71 | #evaluate = True, |
---|
| 72 | dependencies = [filename], |
---|
| 73 | verbose = True) |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | |
---|
| 78 | #Checks |
---|
| 79 | #print reference_time |
---|
| 80 | #print input_time |
---|
| 81 | assert reference_time[0] == 0.0 |
---|
| 82 | assert reference_time[-1] == finaltime |
---|
| 83 | assert allclose(reference_time, input_time) |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | #Validation |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | for k, name in enumerate(gauge_names): |
---|
| 91 | sqsum = 0 |
---|
| 92 | denom = 0 |
---|
| 93 | model = [] |
---|
| 94 | print 'Validating ' + name |
---|
| 95 | for i, t in enumerate(reference_time): |
---|
| 96 | ref = gauge_values[k][i] |
---|
| 97 | val = f(t, point_id = k)[0] |
---|
| 98 | model.append(val) |
---|
| 99 | |
---|
| 100 | sqsum += (ref - val)**2 |
---|
| 101 | denom += ref**2 |
---|
| 102 | |
---|
| 103 | print sqsum |
---|
| 104 | print sqsum/denom |
---|
| 105 | |
---|
| 106 | from pylab import * |
---|
| 107 | ion() |
---|
| 108 | hold(False) |
---|
| 109 | plot(reference_time, gauge_values[k], 'r.-', |
---|
| 110 | reference_time, model, 'k-') |
---|
| 111 | title('Gauge %s' %name) |
---|
| 112 | xlabel('time(s)') |
---|
| 113 | ylabel('stage (m)') |
---|
| 114 | legend(('Observed', 'Modelled'), shadow=True, loc='upper left') |
---|
| 115 | savefig(name, dpi = 300) |
---|
| 116 | |
---|
| 117 | raw_input('Next') |
---|
| 118 | show() |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | #from pylab import * |
---|
| 123 | #plot(time, stage) |
---|
| 124 | #show() |
---|