I am trying to perform a simple 2D linear interpolation with Scipy interpolation.griddata but it behaves in a strange way : it runs forever and the computation can't be interrupted (100% CPU, RAM doesn't move, and I have to kill the process). The data + coordinates that should be interpolated are located inside a polygon. At first, I tried to use a meshgrid containing only points located inside that polygon, but it didn't solve the problem.
I also tried to reduce to a very tiny interpolation grid (10*10), but it kept on freezing.
Then I tried to interpolate the meshgrid points one by one and found out that only some of them were causing the issue, and I can't figure out why.
Here is a sample of the code :
# Sets negative levels to zero
x_a = x_a.clip(min = 0)
coordinate = np.array(coordinate)
# Experiment zone polygon.
if city == "sf":
zone = [[-122.397752, 37.799137],
[-122.423329, 37.795746],
[-122.418523, 37.772886],
[-122.39243, 37.793982]]
zone.append(zone[0]) # to close the polygon.
# Background map extent lat/long coordinates.
x_min_map = min([xy[0] for xy in zone])
x_max_map = max([xy[0] for xy in zone])
y_min_map = min([xy[1] for xy in zone])
y_max_map = max([xy[1] for xy in zone])
def IsPointInsidePolygon(x, y, poly):
n = len(poly)
inside = False
p1x, p1y = poly[0]
for i in range(n + 1):
p2x, p2y = poly[i % n]
if y > min(p1y, p2y):
if y <= max(p1y, p2y):
if x <= max(p1x, p2x):
if p1y != p2y:
xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
if p1x == p2x or x <= xinters:
inside = not inside
p1x, p1y = p2x, p2y
return inside
# Creates a mesh grid adapted to the sf zone, on which to interpolate.
N = 10
xi = np.linspace(x_min_map, x_max_map, N)
yi = np.linspace(y_min_map, y_max_map, N)
xi, yi = np.meshgrid(xi, yi)
xi, yi = xi.flatten(), yi.flatten()
all_points = np.vstack((xi, yi)).T
delete_index = []
for i, point in enumerate(all_points):
if not IsPointInsidePolygon(point[0], point[1], zone):
delete_index.append(i)
grid_points = np.delete(all_points, delete_index, axis = 0)
# Interpolation
zi = interpolate.griddata(coordinate, x_a, grid_points, method = "linear")
The "coordinate" variable is an array that contains approximately 1e4 float coordinates that are located inside the "zone" polygon.
The x_min_map etc. are the computed max and min x, y values of the polygon, so that the original meshgrid is the smallest rectangle containing the whole polygon.
Any advice will be appreciated !
