1 | """Least squares interpolation. |
---|
2 | |
---|
3 | Implements a least-squares interpolation. |
---|
4 | |
---|
5 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
6 | Geoscience Australia, 2004. |
---|
7 | |
---|
8 | DESIGN ISSUES |
---|
9 | * what variables should be global? |
---|
10 | - if there are no global vars functions can be moved around alot easier |
---|
11 | |
---|
12 | * The public interface |
---|
13 | __init__ |
---|
14 | interpolate |
---|
15 | interpolate_block |
---|
16 | |
---|
17 | """ |
---|
18 | |
---|
19 | import time |
---|
20 | import os |
---|
21 | from warnings import warn |
---|
22 | |
---|
23 | from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ |
---|
24 | ArrayType, allclose, take, NewAxis, arange |
---|
25 | |
---|
26 | from anuga.caching.caching import cache |
---|
27 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
28 | from anuga.utilities.sparse import Sparse, Sparse_CSR |
---|
29 | from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError |
---|
30 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
31 | from anuga.utilities.quad import build_quadtree |
---|
32 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
33 | from anuga.utilities.polygon import in_and_outside_polygon |
---|
34 | from anuga.geospatial_data.geospatial_data import Geospatial_data, \ |
---|
35 | ensure_absolute |
---|
36 | from search_functions import search_tree_of_vertices |
---|
37 | |
---|
38 | |
---|
39 | |
---|
40 | class FitInterpolate: |
---|
41 | |
---|
42 | def __init__(self, |
---|
43 | vertex_coordinates, |
---|
44 | triangles, |
---|
45 | mesh_origin=None, |
---|
46 | verbose=False, |
---|
47 | max_vertices_per_cell=30): |
---|
48 | |
---|
49 | |
---|
50 | """ Build interpolation matrix mapping from |
---|
51 | function values at vertices to function values at data points |
---|
52 | |
---|
53 | Inputs: |
---|
54 | |
---|
55 | vertex_coordinates: List of coordinate pairs [xi, eta] of |
---|
56 | points constituting a mesh (or an m x 2 Numeric array or |
---|
57 | a geospatial object) |
---|
58 | Points may appear multiple times |
---|
59 | (e.g. if vertices have discontinuities) |
---|
60 | |
---|
61 | triangles: List of 3-tuples (or a Numeric array) of |
---|
62 | integers representing indices of all vertices in the mesh. |
---|
63 | |
---|
64 | mesh_origin: A geo_reference object or 3-tuples consisting of |
---|
65 | UTM zone, easting and northing. |
---|
66 | If specified vertex coordinates are assumed to be |
---|
67 | relative to their respective origins. |
---|
68 | |
---|
69 | max_vertices_per_cell: Number of vertices in a quad tree cell |
---|
70 | at which the cell is split into 4. |
---|
71 | |
---|
72 | Note: Don't supply a vertex coords as a geospatial object and |
---|
73 | a mesh origin, since geospatial has its own mesh origin. |
---|
74 | """ |
---|
75 | |
---|
76 | #Convert input to Numeric arrays |
---|
77 | triangles = ensure_numeric(triangles, Int) |
---|
78 | vertex_coordinates = ensure_absolute(vertex_coordinates, |
---|
79 | geo_reference = mesh_origin) |
---|
80 | |
---|
81 | #Don't pass geo_reference to mesh. It doesn't work. (Still??) |
---|
82 | |
---|
83 | if verbose: print 'FitInterpolate: Building mesh' |
---|
84 | self.mesh = Mesh(vertex_coordinates, triangles) |
---|
85 | self.mesh.check_integrity() |
---|
86 | |
---|
87 | if verbose: print 'FitInterpolate: Building quad tree' |
---|
88 | self.root = build_quadtree(self.mesh, |
---|
89 | max_points_per_cell = max_vertices_per_cell) |
---|
90 | #print "self.root",self.root.show() |
---|
91 | |
---|
92 | |
---|
93 | def __repr__(self): |
---|
94 | return 'Interpolation object based on: ' + repr(self.mesh) |
---|