I have a 2d array of directional data. I need to interpolate over a higher resolution grid however the ready made functions like scipy interp2d, etc don't account for the discontinuity between 0 and 360.
I have code for doing this for a single grid of 4 points (thanks How to perform bilinear interpolation in Python and Rotation Interpolation) however I would like it to accept big data sets at once - just like the interp2d function. How can I encorporate this into the code below in a way which doesn't just loop over all of the data?
Thanks!
def shortest_angle(beg,end,amount):
shortest_angle=((((end - beg) % 360) + 540) % 360) - 180
return shortest_angle*amount
def bilinear_interpolation_rotation(x, y, points):
'''Interpolate (x,y) from values associated with four points.
The four points are a list of four triplets: (x, y, value).
The four points can be in any order. They should form a rectangle.
'''
points = sorted(points) # order points by x, then by y
(x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
raise ValueError('points do not form a rectangle')
if not x1 <= x <= x2 or not y1 <= y <= y2:
raise ValueError('(x, y) not within the rectangle')
# interpolate over the x value at each y point
fxy1 = q11 + shortest_angle(q11,q21,((x-x1)/(x2-x1)))
fxy2 = q12 + shortest_angle(q12,q22,((x-x1)/(x2-x1)))
# interpolate over the y values
fxy = fxy1 + shortest_angle(fxy1,fxy2,((y-y1)/(y2-y1)))
return fxy
I'm going to reuse some personal Point and Point3D simplified classes for this example:
Pointclass Point:
#Constructors
def __init__(self, x, y):
self.x = x
self.y = y
# Properties
@property
def x(self):
return self._x
@x.setter
def x(self, value):
self._x = float(value)
@property
def y(self):
return self._y
@y.setter
def y(self, value):
self._y = float(value)
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) < (other.x, other.y)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) <= (other.x, other.y)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) > (other.x, other.y)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) >= (other.x, other.y)
It represents a 2D point. It has a simple constructor, x and y properties that ensure they always store floats, magic methods for string representation as (x,y) and comparison to make them sortable (sorts by x, then by y). My original class has additional features such as addition and substraction (vector behaviour) magic methods both but they are not needed for this example.
Point3Dclass Point3D(Point):
# Constructors
def __init__(self, x, y, z):
super().__init__(x, y)
self.z = z
@classmethod
def from2D(cls, p, z):
return cls(p.x, p.y, z)
# Properties
@property
def z(self):
return self._z
@z.setter
def z(self, value):
self._z = (value + 180.0) % 360 - 180
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y},{p.z})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y') and hasattr(other, 'z')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y) and (self.z == other.z)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y) or (self.z != other.z)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) < (other.x, other.y, other.z)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) <= (other.x, other.y, other.z)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) > (other.x, other.y, other.z)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) >= (other.x, other.y, other.z)
Same as Point but for 3D points. It also includes an additional constructor classmethod that takes a Point and its z value as arguments.
def linear_interpolation(x, *points, extrapolate=False):
# Check there are a minimum of two points
if len(points) < 2:
raise ValueError("Not enought points given for interpolation.")
# Sort the points
points = sorted(points)
# Check that x is the valid interpolation interval
if not extrapolate and (x < points[0].x or x > points[-1].x):
raise ValueError("{} is not in the interpolation interval.".format(x))
# Determine which are the two surrounding interpolation points
if x < points[0].x:
i = 0
elif x > points[-1].x:
i = len(points)-2
else:
i = 0
while points[i+1].x < x:
i += 1
p1, p2 = points[i:i+2]
# Interpolate
return Point(x, p1.y + (p2.y-p1.y) * (x-p1.x) / (p2.x-p1.x))
It takes a first position argument that will determine the x whose y value we want to calculate, and an infinite amount of Point instances from where we want to interpolate. A keyword argument (extrapolate) allows to turn on extrapolation. A Point instance is returned with the requested x and the calculated y values.
I offer two alternatives, both of them have a similar signature to the previous interpolation function. A Point whose z value we want to calculate, a keyword argument (extrapolate) that turns on extrapolation and return a Point3D instance with the requested and calculated data. The difference between these two approaches are how the values that will be used to interpolate are provided:
The first approach takes a two-levels-deep nested dict. The first level keys represent the x values, the second level keys the y values and the second level values the z values.
def bilinear_interpolation(p, points, extrapolate=False):
x_values = sorted(points.keys())
# Check there are a minimum of two x values
if len(x_values) < 2:
raise ValueError("Not enought points given for interpolation.")
y_values = set()
for value in points.values():
y_values.update(value.keys())
y_values = sorted(y_values)
# Check there are a minimum of two y values
if len(y_values) < 2:
raise ValueError("Not enought points given for interpolation.")
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for i, surrounding in enumerate(surroundings):
try:
surroundings[i] = Point3D.from2D(surrounding, points[surrounding.x][surrounding.y])
except KeyError:
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), {1: {2: 5, 4: 6}, 3: {2: 3, 4: 9}}))
The second approach takes an infinite number of Point3D instances.
def bilinear_interpolation(p, *points, extrapolate=False):
# Check there are a minimum of four points
if len(points) < 4:
raise ValueError("Not enought points given for interpolation.")
# Sort the points into a grid
x_values = set()
y_values = set()
for point in sorted(points):
x_values.add(point.x)
y_values.add(point.y)
x_values = sorted(x_values)
y_values = sorted(y_values)
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for point in points:
for i, surrounding in enumerate(surroundings):
if point.x == surrounding.x and point.y == surrounding.y:
surroundings[i] = point
for surrounding in surroundings:
if not isinstance(surrounding, Point3D):
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), Point3D(3,2,3), Point3D(1,4,6), Point3D(3,4,9), Point3D(1,2,5)))
You can see from both approaches that they use the previously defined linear_interpoaltion function, and that they always set extrapolation to True as they already raised an exception if it was False and the requested point was outside the provided interval.