2

I built a plot over roads graph and am trying to interpolate the data in a raster, but some nodes cause triangulation error.

Earlier I noticed and fixed a problem when multiple points sticked in one place on a graph and this caused triangulation error. Now these duplicates are removed, but still some sets of points make it crash.

I did not manage to detect what in particular. When I order them by a coordinate and remove either half, the problem disappears.

I suspect some points be standing in one line close enough that the triangulator algorithm tries to build a triangle of them and crashes. But have not found them yet.

What is a way to find the set of points causing the error?

Here's the code, it comes from python-osrm package which copies it from this example in matplotlib:

import csv
import matplotlib
import itertools
import numpy as np
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata


with open('/tmp/triangle.csv') as f:
    rd = csv.reader(f)
    matrix = [map(float, r) for r in rd]

x, y, z = zip(*matrix)

minx = np.nanmin(x)
miny = np.nanmin(y)
maxx = np.nanmax(x)
maxy = np.nanmax(y)

xi = np.linspace(minx, maxx, 400)
yi = np.linspace(miny, maxy, 400)
zi = griddata(x, y, z, xi, yi, interp='linear')

Here's the data: https://pastebin.com/TWz3vNbC

1 Answer 1

2

The solution was to use another interpolation algorithm. CubicTriInterpolator can process a mesh of collinear points.

minx = np.nanmin(x)
miny = np.nanmin(y)
maxx = np.nanmax(x)
maxy = np.nanmax(y)

# Assuming we want a square grid for the interpolation
xi = np.linspace(minx, maxx, 400)
yi = np.linspace(miny, maxy, 400)

# make all combinations of Xi and Yi (400*400), then unwrap it in one line (in 2 rows (x, y) by 160k)
ci = np.reshape(np.meshgrid(xi, yi), (2, 160000))

triang = Triangulation(x, y)
interp = CubicTriInterpolator(triang, z, kind='geom')

# interpolate values, then re-wrap it into 400 by 400 for contourf (it requires a grid like pixels)
zi = np.reshape(interp(*ci), (400, 400))

collec_poly = plt.contourf(
    xi, yi, zi, levels, cmap=plt.cm.rainbow,
    vmax=abs(zi).max(), vmin=-abs(zi).max(), alpha=0.35
    )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.