Discussion:
No subject
Dhananjay
2014-10-17 03:45:20 UTC
Permalink
Hello,

This might be simple, but I guess I am missing something here.
I have data file as follows:

2.1576318858 -1.8651195165 4.2333428278
2.1681875208 -1.9229968780 4.1989176884
2.3387636157 -2.0376253255 2.4460899122
2.1696565965 -2.6186941271 4.4172007912
2.0848862071 -2.1708981985 3.3404520962
2.0824347942 -1.9142798955 3.3629290206
2.0281685821 -1.8103363482 2.5446721669
2.3309993378 -1.8721153619 2.7006893016
2.0957461483 -1.5379071451 4.5228264441
2.2761376261 -2.5935979811 3.9231744717
.
.
.
(total of 200 lines)

Columns 1,2,3 corresponds to x,y,z axis data points.
This is not a continuous data. I wish to make a plot as a 2D with 3rd
dimension (i.e z-axis data) as a color map with color bar on right hand
side.

As a beginner, I tried to follow tutorial with some modification as
follows:
http://matplotlib.org/examples/pylab_examples/tricontour_vs_griddata.html

# Read data from file:
fl1 = open('flooding-psiphi.dat','r').readlines()
xs = ys = zs = []
for line in fl1:
line = line.split()
xs.append(float(line[0]))
ys.append(float(line[1]))
zs.append(float(line[2]))

print xs[0], ys[0], zs[0]
xi = np.mgrid[-5.0:5.0:200j]
yi = np.mgrid[-5.0:5.0:200j]
zi = griddata((x, y), z, (xi, yi), method='cubic')
plt.subplot(221)

plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
plt.contourf(xi, yi, zi, 15, cmap=plt.cm.rainbow,
norm=plt.Normalize(vmax=abs(zi).max(), vmin=-abs(zi).max()))
plt.colorbar() # draw colorbar
plt.plot(x, y, 'ko', ms=3)
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.title('griddata and contour (%d points, %d grid points)' %
(npts, ngridx*ngridy))
#print ('griddata and contour seconds: %f' % (time.clock() - start))
plt.gcf().set_size_inches(6, 6)
plt.show()


However, I failed and getting long error as follows:

QH6154 qhull precision error: initial facet 1 is coplanar with the interior
point
ERRONEOUS FACET:
- f1
- flags: bottom simplicial upperDelaunay flipped
- normal: 0.7071 -0.7071 0
- offset: -0
- vertices: p600(v2) p452(v1) p304(v0)
- neighboring facets: f2 f3 f4

While executing: | qhull d Qz Qbb Qt
Options selected for Qhull 2010.1 2010/01/14:
run-id 1531309415 delaunay Qz-infinity-point Qbbound-last Qtriangulate
_pre-merge _zero-centrum Pgood _max-width 8.8 Error-roundoff 1.2e-14
_one-merge 8.6e-14 _near-inside 4.3e-13 Visible-distance 2.5e-14
U-coplanar-distance 2.5e-14 Width-outside 4.9e-14 _wide-facet 1.5e-13

precision problems (corrected unless 'Q0' or an error)
2 flipped facets

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p228(v3): 2.4 2.4 1.4
- p600(v2): 1.4 1.4 8.8
- p452(v1): 5.7 5.7 8
- p304(v0): -3.1 -3.1 2.4

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet. The maximum round off error for
computing distances is 1.2e-14. The center point, facets and distances
to the center point are as follows:

center point 1.595 1.595 5.173

facet p600 p452 p304 distance= 0
facet p228 p452 p304 distance= 0
facet p228 p600 p304 distance= 0
facet p228 p600 p452 distance= 0

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates. Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
0: -3.134 5.701 difference= 8.835
1: -3.134 5.701 difference= 8.835
2: -2.118e-22 8.835 difference= 8.835

If the input should be full dimensional, you have several options that
may determine an initial simplex:
- use 'QJ' to joggle the input and make it full dimensional
- use 'QbB' to scale the points to the unit cube
- use 'QR0' to randomly rotate the input for different maximum points
- use 'Qs' to search all points for the initial simplex
- use 'En' to specify a maximum roundoff error less than 1.2e-14.
- trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
- use 'QJ' to joggle the input and make it full dimensional
- use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
pick the coordinate with the least range. The hull will have the
correct topology.
- determine the flat containing the points, rotate the points
into a coordinate plane, and delete the other coordinates.
- add one or more points to make the input full dimensional.
Traceback (most recent call last):
File "./scatter.py", line 43, in <module>
zi = griddata((x, y), z, (xi, yi), method='linear')
File "/usr/lib/python2.7/dist-packages/scipy/interpolate/ndgriddata.py",
line 183, in griddata
ip = LinearNDInterpolator(points, values, fill_value=fill_value)
File "interpnd.pyx", line 192, in
scipy.interpolate.interpnd.LinearNDInterpolator.__init__
(scipy/interpolate/interpnd.c:2598)
File "qhull.pyx", line 948, in scipy.spatial.qhull.Delaunay.__init__
(scipy/spatial/qhull.c:4121)
File "qhull.pyx", line 172, in scipy.spatial.qhull._construct_delaunay
(scipy/spatial/qhull.c:1314)
RuntimeError: Qhull error


Could anyone help me to point out what exactly I am missing here.
I just wish to plot 2D map with color bar for Z-axis.

Thank you in advance.


-- DJ
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-list/attachments/20141017/e95c9fe6/attachment.html>
Cameron Simpson
2014-10-17 10:43:16 UTC
Permalink
Post by Dhananjay
This might be simple, but I guess I am missing something here.
2.1576318858 -1.8651195165 4.2333428278
2.1681875208 -1.9229968780 4.1989176884
2.3387636157 -2.0376253255 2.4460899122
2.1696565965 -2.6186941271 4.4172007912
2.0848862071 -2.1708981985 3.3404520962
2.0824347942 -1.9142798955 3.3629290206
2.0281685821 -1.8103363482 2.5446721669
2.3309993378 -1.8721153619 2.7006893016
2.0957461483 -1.5379071451 4.5228264441
2.2761376261 -2.5935979811 3.9231744717
(total of 200 lines)
Columns 1,2,3 corresponds to x,y,z axis data points.
This is not a continuous data. I wish to make a plot as a 2D with 3rd
dimension (i.e z-axis data) as a color map with color bar on right hand side.
http://matplotlib.org/examples/pylab_examples/tricontour_vs_griddata.html
[...]

Initially I've just got comments about the code in general, though they may
help you find your issues.
Post by Dhananjay
fl1 = open('flooding-psiphi.dat','r').readlines()
This reads all the lines into memory, into the list "fl1".
For 200 lines this is not a bit issue, but since you process them immediately
you are usually better off reading the file progressively.
Post by Dhananjay
xs = ys = zs = []
This is actually a bug in your code. I would guess you've based it on other
example code such as:

x = y = z = 0

The problem is that both these assignment statements assign the _same_ object
to all 3 variables. With an int or a str this isn't an issue because they are
immutable; any further math or modification actually makes new objects.

However, you use:

xs = ys = zs = []

This makes exactly _one_ list, and assigns it ("binds it" in Python terms) to
all three names.

The result is that when you append to the list, you're appending to the same
list every time. Effectively this means (1) that xs, ys and zs have the same
values and (2) they're 3 times as long as they should be.

Do this:

xs = []
ys = []
zs = []

That makes three separate lists.
Post by Dhananjay
??? line = line.split()
??? xs.append(float(line[0]))
??? ys.append(float(line[1]))
??? zs.append(float(line[2]))
Returning to the "reading the file progressively" thing, for a really big file
(and as a matter of general practice) you're better off writing this like:

for line in open('flooding-psiphi.dat','r'):
??? line = line.split()
??? xs.append(float(line[0]))
??? ys.append(float(line[1]))
??? zs.append(float(line[2]))

which only ever reads one line of text from the file at a time.

Start by fixing the:

xs = ys = zs = []

assignment and see how the behaviour changes.

Cheers,
Cameron Simpson <cs at zip.com.au>
Terry Reedy
2014-10-17 16:37:03 UTC
Permalink
Post by Dhananjay
Post by Dhananjay
2.1576318858 -1.8651195165 4.2333428278
...
(total of 200 lines)
Columns 1,2,3 corresponds to x,y,z axis data points.
line = line.split()
xs.append(float(line[0]))
ys.append(float(line[1]))
zs.append(float(line[2]))
A further refinement:
for line in open('flooding-psiphi.dat','r'):
x, y, z = map(float, line.split())
xs.append(x)
ys.append(y)
zs.append(z)
--
Terry Jan Reedy
Continue reading on narkive:
Loading...