| 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 anuga.caching.caching import cache |
|---|
| 24 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
|---|
| 25 | |
|---|
| 26 | from anuga.utilities.sparse import Sparse, Sparse_CSR |
|---|
| 27 | from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError |
|---|
| 28 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 29 | from anuga.utilities.polygon import in_and_outside_polygon |
|---|
| 30 | from anuga.utilities.quad import build_quadtree |
|---|
| 31 | |
|---|
| 32 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
|---|
| 33 | from anuga.geospatial_data.geospatial_data import Geospatial_data, \ |
|---|
| 34 | ensure_absolute |
|---|
| 35 | from anuga.fit_interpolate.search_functions import set_last_triangle |
|---|
| 36 | |
|---|
| 37 | import numpy as num |
|---|
| 38 | |
|---|
| 39 | |
|---|
| 40 | # tests fail if 2 is used |
|---|
| 41 | MAX_VERTICES_PER_CELL = 13 # A value of 8 or lower can cause problems, |
|---|
| 42 | # if a vert has 9 triangles. |
|---|
| 43 | |
|---|
| 44 | build_quadtree_time = 0 |
|---|
| 45 | |
|---|
| 46 | def get_build_quadtree_time(): |
|---|
| 47 | return build_quadtree_time |
|---|
| 48 | |
|---|
| 49 | class FitInterpolate: |
|---|
| 50 | |
|---|
| 51 | def __init__(self, |
|---|
| 52 | vertex_coordinates=None, |
|---|
| 53 | triangles=None, |
|---|
| 54 | mesh=None, |
|---|
| 55 | mesh_origin=None, |
|---|
| 56 | verbose=False, |
|---|
| 57 | max_vertices_per_cell=None): |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | """ Build interpolation matrix mapping from |
|---|
| 61 | function values at vertices to function values at data points |
|---|
| 62 | |
|---|
| 63 | Pass in a mesh instance or vertex_coordinates and triangles |
|---|
| 64 | and optionally mesh_origin |
|---|
| 65 | |
|---|
| 66 | Inputs: |
|---|
| 67 | |
|---|
| 68 | vertex_coordinates: List of coordinate pairs [xi, eta] of |
|---|
| 69 | points constituting a mesh (or an m x 2 numeric array or |
|---|
| 70 | a geospatial object) |
|---|
| 71 | Points may appear multiple times |
|---|
| 72 | (e.g. if vertices have discontinuities) |
|---|
| 73 | |
|---|
| 74 | triangles: List of 3-tuples (or a numeric array) of |
|---|
| 75 | integers representing indices of all vertices in the mesh. |
|---|
| 76 | |
|---|
| 77 | mesh: A mesh instance describing the mesh. |
|---|
| 78 | |
|---|
| 79 | mesh_origin: A geo_reference object or 3-tuples consisting of |
|---|
| 80 | UTM zone, easting and northing. |
|---|
| 81 | If specified vertex coordinates are assumed to be |
|---|
| 82 | relative to their respective origins. |
|---|
| 83 | |
|---|
| 84 | max_vertices_per_cell: Number of vertices in a quad tree cell |
|---|
| 85 | at which the cell is split into 4. |
|---|
| 86 | |
|---|
| 87 | Note: Don't supply a vertex coords as a geospatial object and |
|---|
| 88 | a mesh origin, since geospatial has its own mesh origin. |
|---|
| 89 | """ |
|---|
| 90 | global build_quadtree_time |
|---|
| 91 | if max_vertices_per_cell == None: |
|---|
| 92 | max_vertices_per_cell = MAX_VERTICES_PER_CELL |
|---|
| 93 | if mesh is None: |
|---|
| 94 | if vertex_coordinates is not None and triangles is not None: |
|---|
| 95 | # Fixme (DSG) Throw errors if triangles or vertex_coordinates |
|---|
| 96 | # are None |
|---|
| 97 | |
|---|
| 98 | #Convert input to numeric arrays |
|---|
| 99 | triangles = ensure_numeric(triangles, num.int) |
|---|
| 100 | vertex_coordinates = ensure_absolute(vertex_coordinates, |
|---|
| 101 | geo_reference = mesh_origin) |
|---|
| 102 | |
|---|
| 103 | if verbose: print 'FitInterpolate: Building mesh' |
|---|
| 104 | self.mesh = Mesh(vertex_coordinates, triangles) |
|---|
| 105 | #self.mesh.check_integrity() # Time consuming |
|---|
| 106 | else: |
|---|
| 107 | self.mesh = None |
|---|
| 108 | else: |
|---|
| 109 | self.mesh = mesh |
|---|
| 110 | |
|---|
| 111 | if self.mesh is not None: |
|---|
| 112 | if verbose: print 'FitInterpolate: Building quad tree' |
|---|
| 113 | #This stores indices of vertices |
|---|
| 114 | t0 = time.time() |
|---|
| 115 | #print "self.mesh.get_extent(absolute=True)", \ |
|---|
| 116 | #self.mesh.get_extent(absolute=True) |
|---|
| 117 | self.root = build_quadtree(self.mesh, |
|---|
| 118 | max_points_per_cell = max_vertices_per_cell) |
|---|
| 119 | #print "self.root",self.root.show() |
|---|
| 120 | |
|---|
| 121 | build_quadtree_time = time.time()-t0 |
|---|
| 122 | set_last_triangle() |
|---|
| 123 | |
|---|
| 124 | def __repr__(self): |
|---|
| 125 | return 'Interpolation object based on: ' + repr(self.mesh) |
|---|
