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,argmin |
---|
10 | from pylab import plot, ion, hold,savefig,semilogx,plotting,loglog |
---|
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=1000, |
---|
21 | filename=mesh_dir_name, |
---|
22 | use_cache=True, |
---|
23 | verbose=True) |
---|
24 | |
---|
25 | |
---|
26 | #topo_dir_name = 'pt_hedland_export.txt' |
---|
27 | topo_dir_name = 'pt_hedland_small.txt' |
---|
28 | #sample = 'pt_hedland_small.txt' |
---|
29 | remainder ='remainder.txt' |
---|
30 | #split topo data |
---|
31 | G = Geospatial_data(file_name = topo_dir_name) |
---|
32 | print 'start split' |
---|
33 | |
---|
34 | |
---|
35 | |
---|
36 | |
---|
37 | |
---|
38 | G_sample,G_loss= G.split(0.01, True) |
---|
39 | |
---|
40 | #G_sample.export_points_file(sample) |
---|
41 | |
---|
42 | G_small, G_other = G_sample.split(0.1,True) |
---|
43 | |
---|
44 | #import sys |
---|
45 | #sys.exit() |
---|
46 | |
---|
47 | alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0,1000.0,10000.0] |
---|
48 | #alphas = [0.001,1,100] |
---|
49 | domains = {} |
---|
50 | #domains = [] |
---|
51 | |
---|
52 | print 'Setup computational domain' |
---|
53 | for alpha in alphas: |
---|
54 | |
---|
55 | #think about making domain absolute when it is created |
---|
56 | domain = Domain(mesh_dir_name, use_cache=False, verbose=True) |
---|
57 | print domain.statistics() |
---|
58 | domain.set_quantity('elevation', |
---|
59 | geospatial_data = G_other, |
---|
60 | use_cache = True, |
---|
61 | verbose = True, |
---|
62 | alpha = alpha) |
---|
63 | domains[alpha]=domain |
---|
64 | |
---|
65 | points=G_small.get_data_points() |
---|
66 | data=array([],typecode=Float) |
---|
67 | data=resize(data,(len(points),3+len(alphas))) |
---|
68 | #gets relative point from sample |
---|
69 | |
---|
70 | data[:,0]=points[:,0] |
---|
71 | data[:,1]=points[:,1] |
---|
72 | print'geo',domain.geo_reference,points[0:10],points[0:10,0] |
---|
73 | elevation_sample=G_small.get_attributes(attribute_name='elevation') |
---|
74 | |
---|
75 | data[:,2]=elevation_sample |
---|
76 | |
---|
77 | print'data',data[0:10],len(alphas) |
---|
78 | normal_cov=array(zeros([len(alphas),2]),typecode=Float) |
---|
79 | i=0 |
---|
80 | |
---|
81 | for domain in domains: |
---|
82 | # print'domain',domain |
---|
83 | points_geo=domains[domain].geo_reference.change_points_geo_ref(points) |
---|
84 | |
---|
85 | |
---|
86 | print'poin',points_geo[0:10],elevation_sample[0:10] |
---|
87 | print'quantities',domains[domain].get_quantity_names() |
---|
88 | elevation_predicted=domains[domain].quantities['elevation'].get_values(interpolation_points=points_geo) |
---|
89 | print'elevation_predicted',elevation_predicted[0:10],len(points),len(elevation_predicted) |
---|
90 | |
---|
91 | data[:,i+3]=elevation_predicted |
---|
92 | sample_cov= cov(elevation_sample) |
---|
93 | |
---|
94 | ele_cov= cov(elevation_sample-elevation_predicted) |
---|
95 | normal_cov[i,:]= [domain,ele_cov/sample_cov] |
---|
96 | |
---|
97 | print'cov',alpha,'= ',normal_cov, ele_cov, sample_cov |
---|
98 | i+=1 |
---|
99 | |
---|
100 | |
---|
101 | print'data',data[100:200] |
---|
102 | print'normal cov',normal_cov |
---|
103 | #to sort array by column |
---|
104 | normal_cov0=normal_cov[:,0] |
---|
105 | normal_cov_new=take(normal_cov,argsort(normal_cov0)) |
---|
106 | semilogx(normal_cov_new[:,0],normal_cov_new[:,1]) |
---|
107 | #loglog(normal_cov_new[:,0],normal_cov_new[:,1]) |
---|
108 | print 'normal_cov_new',normal_cov_new |
---|
109 | savefig("alphas",dpi=300) |
---|
110 | |
---|
111 | print argmin(normal_cov_new,axis=1) |
---|
112 | print min(normal_cov_new[:,1]) , normal_cov_new[(argmin(normal_cov_new,axis=0))[1],0] |
---|
113 | |
---|