1 | |
---|
2 | from anuga.abstract_2d_finite_volumes import domain |
---|
3 | from anuga.shallow_water import Domain |
---|
4 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
5 | from anuga.coordinate_transforms import Geo_reference |
---|
6 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
7 | from anuga.shallow_water import Reflective_boundary |
---|
8 | from anuga.utilities.numerical_tools import cov |
---|
9 | from Numeric import array, resize,shape,Float,zeros,take,argsort |
---|
10 | from pylab import plot, ion, hold,savefig,semilogx,plotting |
---|
11 | |
---|
12 | |
---|
13 | poly_topo = [[664000,7752000],[664000,7754000], |
---|
14 | [666000,7754000],[666000,7752000]] |
---|
15 | mesh_dir_name = 'temp.msh' |
---|
16 | |
---|
17 | create_mesh_from_regions(poly_topo, |
---|
18 | boundary_tags={'back': [2], 'side': [1,3], |
---|
19 | 'ocean': [0]}, |
---|
20 | maximum_triangle_area=100, |
---|
21 | filename=mesh_dir_name, |
---|
22 | use_cache=True, |
---|
23 | verbose=True) |
---|
24 | |
---|
25 | |
---|
26 | topo_dir_name = 'pt_hedland_export.txt' |
---|
27 | sample = 'sample.txt' |
---|
28 | remainder ='remainder.txt' |
---|
29 | #split topo data |
---|
30 | G = Geospatial_data(file_name = topo_dir_name) |
---|
31 | print 'start split' |
---|
32 | |
---|
33 | G_sample,G_loss= G.split(0.1, True) |
---|
34 | |
---|
35 | G_small, G_other = G_sample.split(0.1,True) |
---|
36 | |
---|
37 | |
---|
38 | |
---|
39 | alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0,1000.0,10000.0] |
---|
40 | #alphas = [0.001,1,100] |
---|
41 | domains = {} |
---|
42 | #domains = [] |
---|
43 | |
---|
44 | print 'Setup computational domain' |
---|
45 | for alpha in alphas: |
---|
46 | |
---|
47 | #think about making domain absolute when it is created |
---|
48 | domain = Domain(mesh_dir_name, use_cache=False, verbose=True) |
---|
49 | print domain.statistics() |
---|
50 | domain.set_quantity('elevation', |
---|
51 | geospatial_data = G_other, |
---|
52 | use_cache = True, |
---|
53 | verbose = True, |
---|
54 | alpha = alpha) |
---|
55 | domains[alpha]=domain |
---|
56 | |
---|
57 | points=G_small.get_data_points() |
---|
58 | data=array([],typecode=Float) |
---|
59 | data=resize(data,(len(points),3+len(alphas))) |
---|
60 | #gets relative point from sample |
---|
61 | |
---|
62 | data[:,0]=points[:,0] |
---|
63 | data[:,1]=points[:,1] |
---|
64 | print'geo',domain.geo_reference,points[0:10],points[0:10,0] |
---|
65 | elevation_sample=G_small.get_attributes(attribute_name='elevation') |
---|
66 | |
---|
67 | data[:,2]=elevation_sample |
---|
68 | |
---|
69 | print'data',data[0:10],len(alphas) |
---|
70 | normal_cov=array(zeros([len(alphas),2]),typecode=Float) |
---|
71 | i=0 |
---|
72 | |
---|
73 | for domain in domains: |
---|
74 | # print'domain',domain |
---|
75 | points_geo=domains[domain].geo_reference.change_points_geo_ref(points) |
---|
76 | |
---|
77 | |
---|
78 | print'poin',points_geo[0:10],elevation_sample[0:10] |
---|
79 | print'quantities',domains[domain].get_quantity_names() |
---|
80 | elevation_predicted=domains[domain].quantities['elevation'].get_values(interpolation_points=points_geo) |
---|
81 | print'elevation_predicted',elevation_predicted[0:10],len(points),len(elevation_predicted) |
---|
82 | |
---|
83 | data[:,i+3]=elevation_predicted |
---|
84 | sample_cov= cov(elevation_sample) |
---|
85 | |
---|
86 | ele_cov= cov(elevation_sample-elevation_predicted) |
---|
87 | normal_cov[i,:]= [domain,ele_cov/sample_cov] |
---|
88 | |
---|
89 | print'cov',alpha,'= ',normal_cov, ele_cov, sample_cov |
---|
90 | i+=1 |
---|
91 | |
---|
92 | |
---|
93 | print'data',data[100:200] |
---|
94 | print'normal cov',normal_cov |
---|
95 | normal_cov0=normal_cov[:,0] |
---|
96 | normal_cov_new=take(normal_cov,argsort(normal_cov0)) |
---|
97 | semilogx(normal_cov_new[:,0],normal_cov_new[:,1]) |
---|
98 | savefig("alphas",dpi=300) |
---|
99 | |
---|