Interpolation is added.

chrest/ablateData.py

@@ -8,6 +8,7 @@
 from chrestData import ChrestData
 from supportPaths import expand_path
 from scipy.spatial import KDTree
+from interpolTree import *
 
 
 class AblateData:
@@ -218,16 +219,21 @@ def map_to_chrest_data(self, chrest_data, field_mapping, field_select_components
         # now search and copy over data
         tree = KDTree(cell_centers)
         dist, points = tree.query(chrest_cell_centers)
+        dist1, points_ablate = tree.query(cell_centers)
 
         # march over each field
         for f in range(len(ablate_field_data)):
             # get in the correct order
-            ablate_field_in_chrest_order = ablate_field_data[f][:, points]
+            ablate_field_in_chrest_order = ablate_field_data[f][:, points_ablate]
 
             setToZero = np.where(dist > max_distance)
 
             ablate_field_in_chrest_order[:, setToZero] = 0.0
 
+            interpolTree = InterpolTree(cell_centers, ablate_field_in_chrest_order.reshape(
+                len(cell_centers), 2 * f + 1))
+            ablate_field_in_chrest_order = interpolTree(chrest_cell_centers, n_near=args.n_near)
+
             # reshape it back to k,j,i
             ablate_field_in_chrest_order = ablate_field_in_chrest_order.reshape(
                 chrest_field_data[f].shape
@@ -274,8 +280,12 @@ def map_to_chrest_data(self, chrest_data, field_mapping, field_select_components
                         )
 
     parser.add_argument('--max_distance', dest='max_distance', type=float,
-                        help="The max distance to search for a point in ablate", default=sys.float_info.max)
+                        help="The max distance to search for a point in ablate", default=sys.float_info.max
+                        )
 
+    parser.add_argument('--n_near', dest='n_near', type=int,
+                        help="number of neighbour cells for interpolation", default=1
+                        )
     args = parser.parse_args()
 
     # this is some example code for chest file post-processing

chrest/interpolTree.py

@@ -0,0 +1,49 @@
+"""
+inverse-distance-weighted interpolation using KDTree:
+"""
+from __future__ import division
+import numpy as np
+from scipy.spatial import cKDTree as KDTree
+
+
+# http://docs.scipy.org/doc/scipy/reference/spatial.html
+# ...............................................................................
+class InterpolTree:
+
+    def __init__(self, X, Y, leafsize=10, stat=0):
+        assert len(X) == len(Y), "len(X) %d != len(z) %d" % (len(X), len(Y))
+        self.tree = KDTree(X, leafsize=leafsize)  # build the tree
+        self.Y = Y
+        self.stat = stat
+        self.wn = 0
+        self.wsum = None;
+
+    def __call__(self, q, n_near=6, eps=.1, p=3, weights=None):
+        # nnear nearest neighbours of each q point --
+        q = np.asarray(q)
+        qdim = q.ndim
+        if qdim == 1:
+            q = np.array([q])
+        if self.wsum is None:
+            self.wsum = np.zeros(n_near)
+
+        self.distances, self.ix = self.tree.query(q, k=n_near, eps=eps)
+        interpol = np.zeros((len(self.distances),) + np.shape(self.Y[0]))
+        jinterpol = 0
+        for dist, ix in zip(self.distances, self.ix):
+            if n_near == 1:
+                wY = self.Y[ix]
+            elif dist[0] < 1e-10:
+                wY = self.Y[ix[0]]
+            else:  # weight Y s by 1/dist --
+                w = 1 / dist ** p
+                if weights is not None:
+                    w *= weights[ix]  # >= 0
+                w /= np.sum(w)
+                wY = np.dot(w, self.Y[ix])
+                if self.stat:
+                    self.wn += 1
+                    self.wsum += w
+            interpol[jinterpol] = wY
+            jinterpol += 1
+        return interpol if qdim > 1 else interpol[0]