|
| 1 | +""" |
| 2 | +Contour plots of unstructured triangular grids. |
| 3 | +""" |
| 4 | +import matplotlib.pyplot as plt |
| 5 | +import matplotlib.tri as tri |
| 6 | +import numpy as np |
| 7 | +import math |
| 8 | + |
| 9 | +# Creating a Triangulation without specifying the triangles results in the |
| 10 | +# Delaunay triangulation of the points. |
| 11 | + |
| 12 | +# First create the x and y coordinates of the points. |
| 13 | +n_angles = 48 |
| 14 | +n_radii = 8 |
| 15 | +min_radius = 0.25 |
| 16 | +radii = np.linspace(min_radius, 0.95, n_radii) |
| 17 | + |
| 18 | +angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False) |
| 19 | +angles = np.repeat(angles[...,np.newaxis], n_radii, axis=1) |
| 20 | +angles[:,1::2] += math.pi/n_angles |
| 21 | + |
| 22 | +x = (radii*np.cos(angles)).flatten() |
| 23 | +y = (radii*np.sin(angles)).flatten() |
| 24 | +z = (np.cos(radii)*np.cos(angles*3.0)).flatten() |
| 25 | + |
| 26 | +# Create the Triangulation; no triangles so Delaunay triangulation created. |
| 27 | +triang = tri.Triangulation(x, y) |
| 28 | + |
| 29 | +# Mask off unwanted triangles. |
| 30 | +xmid = x[triang.triangles].mean(axis=1) |
| 31 | +ymid = y[triang.triangles].mean(axis=1) |
| 32 | +mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0) |
| 33 | +triang.set_mask(mask) |
| 34 | + |
| 35 | +# pcolor plot. |
| 36 | +plt.figure() |
| 37 | +plt.gca().set_aspect('equal') |
| 38 | +plt.tricontourf(triang, z) |
| 39 | +plt.colorbar() |
| 40 | +plt.tricontour(triang, z, colors='k') |
| 41 | +plt.title('Contour plot of Delaunay triangulation') |
| 42 | + |
| 43 | +# You can specify your own triangulation rather than perform a Delaunay |
| 44 | +# triangulation of the points, where each triangle is given by the indices of |
| 45 | +# the three points that make up the triangle, ordered in either a clockwise or |
| 46 | +# anticlockwise manner. |
| 47 | + |
| 48 | +xy = np.asarray([ |
| 49 | + [-0.101,0.872],[-0.080,0.883],[-0.069,0.888],[-0.054,0.890],[-0.045,0.897], |
| 50 | + [-0.057,0.895],[-0.073,0.900],[-0.087,0.898],[-0.090,0.904],[-0.069,0.907], |
| 51 | + [-0.069,0.921],[-0.080,0.919],[-0.073,0.928],[-0.052,0.930],[-0.048,0.942], |
| 52 | + [-0.062,0.949],[-0.054,0.958],[-0.069,0.954],[-0.087,0.952],[-0.087,0.959], |
| 53 | + [-0.080,0.966],[-0.085,0.973],[-0.087,0.965],[-0.097,0.965],[-0.097,0.975], |
| 54 | + [-0.092,0.984],[-0.101,0.980],[-0.108,0.980],[-0.104,0.987],[-0.102,0.993], |
| 55 | + [-0.115,1.001],[-0.099,0.996],[-0.101,1.007],[-0.090,1.010],[-0.087,1.021], |
| 56 | + [-0.069,1.021],[-0.052,1.022],[-0.052,1.017],[-0.069,1.010],[-0.064,1.005], |
| 57 | + [-0.048,1.005],[-0.031,1.005],[-0.031,0.996],[-0.040,0.987],[-0.045,0.980], |
| 58 | + [-0.052,0.975],[-0.040,0.973],[-0.026,0.968],[-0.020,0.954],[-0.006,0.947], |
| 59 | + [ 0.003,0.935],[ 0.006,0.926],[ 0.005,0.921],[ 0.022,0.923],[ 0.033,0.912], |
| 60 | + [ 0.029,0.905],[ 0.017,0.900],[ 0.012,0.895],[ 0.027,0.893],[ 0.019,0.886], |
| 61 | + [ 0.001,0.883],[-0.012,0.884],[-0.029,0.883],[-0.038,0.879],[-0.057,0.881], |
| 62 | + [-0.062,0.876],[-0.078,0.876],[-0.087,0.872],[-0.030,0.907],[-0.007,0.905], |
| 63 | + [-0.057,0.916],[-0.025,0.933],[-0.077,0.990],[-0.059,0.993] ]) |
| 64 | +x = xy[:,0]*180/3.14159 |
| 65 | +y = xy[:,1]*180/3.14159 |
| 66 | +x0 = -5 |
| 67 | +y0 = 52 |
| 68 | +z = np.exp(-0.01*( (x-x0)*(x-x0) + (y-y0)*(y-y0) )) |
| 69 | + |
| 70 | +triangles = np.asarray([ |
| 71 | + [67,66, 1],[65, 2,66],[ 1,66, 2],[64, 2,65],[63, 3,64],[60,59,57], |
| 72 | + [ 2,64, 3],[ 3,63, 4],[ 0,67, 1],[62, 4,63],[57,59,56],[59,58,56], |
| 73 | + [61,60,69],[57,69,60],[ 4,62,68],[ 6, 5, 9],[61,68,62],[69,68,61], |
| 74 | + [ 9, 5,70],[ 6, 8, 7],[ 4,70, 5],[ 8, 6, 9],[56,69,57],[69,56,52], |
| 75 | + [70,10, 9],[54,53,55],[56,55,53],[68,70, 4],[52,56,53],[11,10,12], |
| 76 | + [69,71,68],[68,13,70],[10,70,13],[51,50,52],[13,68,71],[52,71,69], |
| 77 | + [12,10,13],[71,52,50],[71,14,13],[50,49,71],[49,48,71],[14,16,15], |
| 78 | + [14,71,48],[17,19,18],[17,20,19],[48,16,14],[48,47,16],[47,46,16], |
| 79 | + [16,46,45],[23,22,24],[21,24,22],[17,16,45],[20,17,45],[21,25,24], |
| 80 | + [27,26,28],[20,72,21],[25,21,72],[45,72,20],[25,28,26],[44,73,45], |
| 81 | + [72,45,73],[28,25,29],[29,25,31],[43,73,44],[73,43,40],[72,73,39], |
| 82 | + [72,31,25],[42,40,43],[31,30,29],[39,73,40],[42,41,40],[72,33,31], |
| 83 | + [32,31,33],[39,38,72],[33,72,38],[33,38,34],[37,35,38],[34,38,35], |
| 84 | + [35,37,36] ]) |
| 85 | + |
| 86 | +# Rather than create a Triangulation object, can simply pass x, y and triangles |
| 87 | +# arrays to tripcolor directly. It would be better to use a Triangulation object |
| 88 | +# if the same triangulation was to be used more than once to save duplicated |
| 89 | +# calculations. |
| 90 | +plt.figure() |
| 91 | +plt.gca().set_aspect('equal') |
| 92 | +plt.tricontourf(x, y, triangles, z) |
| 93 | +plt.colorbar() |
| 94 | +plt.title('Contour plot of user-specified triangulation') |
| 95 | +plt.xlabel('Longitude (degrees)') |
| 96 | +plt.ylabel('Latitude (degrees)') |
| 97 | + |
| 98 | +plt.show() |
0 commit comments