From c25e5fba7094d79522e92760730211514b2a84f1 Mon Sep 17 00:00:00 2001 From: Dan Sandiford Date: Sun, 7 May 2017 23:47:52 +1000 Subject: [PATCH] nearest neighbour / IDW interpolation using scipy kdTree --- unsupported/__init__.py | 23 ---------------- unsupported/utilities/__init__.py | 13 +++++++++ unsupported/utilities/interpolation.py | 38 ++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 23 deletions(-) delete mode 100644 unsupported/__init__.py create mode 100644 unsupported/utilities/__init__.py create mode 100644 unsupported/utilities/interpolation.py diff --git a/unsupported/__init__.py b/unsupported/__init__.py deleted file mode 100644 index 8700eca04..000000000 --- a/unsupported/__init__.py +++ /dev/null @@ -1,23 +0,0 @@ -##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~## -## ## -## This file forms part of the Underworld geophysics modelling application. ## -## ## -## For full license and copyright information, please refer to the LICENSE.md file ## -## located at the project root, or contact the authors. ## -## ## -##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~## - -""" -Modules found in 'unsupported' have been contributed to the underworld2 code. -These modules are not covered by the testing system, may contain additional python -dependencies, are not yet well documented, are not guaranteed to work with every -mesh type, may be limited to special cases (e.g. only 2D, only Cartesian) and -may never become supported. - - - interfaces: tracking of interfaces defined by point collections -""" - -__version__ = "0.0" - - -# import interfaces diff --git a/unsupported/utilities/__init__.py b/unsupported/utilities/__init__.py new file mode 100644 index 000000000..56ae55988 --- /dev/null +++ b/unsupported/utilities/__init__.py @@ -0,0 +1,13 @@ +__version__ = "0.1" + + +import warnings + +warnings.warn( +"""\n +The alchemy module is not supported. +Questions should be addressed to sandd@student.unimelb.edu.au \n """ +) + + +from interpolation import * diff --git a/unsupported/utilities/interpolation.py b/unsupported/utilities/interpolation.py new file mode 100644 index 000000000..b8b13446d --- /dev/null +++ b/unsupported/utilities/interpolation.py @@ -0,0 +1,38 @@ +from scipy.spatial import cKDTree as kdTree +import numpy as np + +def nn_evaluation(_fromCoords, _toCoords, n=1, weighted=False): + + """ + This function provides nearest neighbour information for uw swarms, + given the "_toCoords", which could be the .data handle (coordinates) of a mesh or a different swarm, + this function returns the indices of the n nearest neighbours in "_fromCoords" (will usually be swarm.particleCoordinates.data ) + it also returns the inverse-distance weights if weighted=True. + + The function works in parallel, if the example below is followed + + + Usage + ------------ + #get the n indexes, weights and distances + ix, weights, d = nn_evaluation(swarm.particleCoordinates.data, toSwarm.particleCoordinates.data, n=n, weighted=False) + + #apply to the 'toSwarm' variable in a parallel-safe way + if len(weights): #parallel safety + toSwarmVar.data[:,0] = np.average(fromSwarmVar.evaluate(fromSwarm)[:,0][ix], weights=weights, axis=len((weights.shape)) - 1) + + """ + + if len(_toCoords) > 0: #this is required for safety in parallel + + tree = kdTree(_fromCoords) + d, ix = tree.query(_toCoords, n) + if n == 1: + weights = np.ones(_toCoords.shape[0]) + elif not weighted: + weights = np.ones((_toCoords.shape[0], n))*(1./n) + else: + weights = (1./d[:])/(1./d[:]).sum(axis=1)[:,None] + return ix, weights, d + else: + return np.empty(0, dtype="int"), np.empty(0, dtype="int"), np.empty(0, dtype="int")