What is the best way to get a visual appealing plot of a singular vector field (if you want to visualize also the field strength). As an example I am playing with the electric fields of two point charges as in the following example: from pylab import *from scipy.integrate import odeintfrom matplotlib import animationfrom matplotlib import cmimport numpy as nprc('font', **{'family': 'serif', 'serif': ['Computer Modern']})rc('text', usetex=True)## Set up chargesclass charge:    def __init__(self, q, pos):        self.q=q        self.pos=poschargesPlus=[]chargesMinus=[]#for i in arange(0,1,1):chargesPlus.append(charge(1,[3,0]))chargesMinus.append(charge(-1,[-3,0]))charges = chargesPlus + chargesMinusdef E_point_charge(q, a, x, y,r):    return q*(x-a[0])/((x-a[0])**2+(y-a[1])**2)**(1.5), \        q*(y-a[1])/((x-a[0])**2+(y-a[1])**2)**(1.5)def E_total(x, y, charges):    Ex, Ey=0, 0    for C in charges:        E=E_point_charge(C.q, C.pos, x, y,1)        Ex=Ex+E[0]        Ey=Ey+E[1]    return [Ex, Ey]domain =2## Cut Quiver plotdef cut(r):    if r < domain:        return 0    else:        return 1cutv = np.vectorize(cut)def cut_total(charges,x):     c = 1    for C in charges:         r = sqrt((C.pos[0] - x[0])**2 + (C.pos[1] - x[1])**2)        c = c*cutv(r)        print c        print C.pos[0],C.pos[1]    return cfig = figure()ax = fig.add_subplot(1,1,1)xMin,xMax=-15,15yMin,yMax=-10,10#ax.plot(x,y)ax.axis('tight')xlim([xMin,xMax])ylim([yMin,yMax])# plot point chargesfor C in charges:    if C.q>0:        plot(C.pos[0], C.pos[1], 'bo', ms=8*sqrt(C.q))    if C.q<0:        plot(C.pos[0], C.pos[1], 'ro', ms=8*sqrt(-C.q))xG,yG = meshgrid(linspace(xMin,xMax,25),linspace(yMin,yMax,25))# plot vector fieldE_totalX,E_totalY = E_total(xG,yG,charges)EAbs = (E_totalX**2 + E_totalY**2)**(0.5)E_XX = E_totalX/EAbsE_YY = E_totalY/EAbs#EAbs = np.nan_to_num(EAbs) #ax.streamplot(xG,yG,E_XX,E_YY,color=EAbs,cmap=cm.autumn)ax.quiver(xG,yG,E_XX,E_YY,EAbs,cmap=cm.GnBu)xlabel('x')ylabel('y')ax.set_aspect(1)plt.savefig('fig1.png')E_totalX = E_totalX*cut_total(charges,[xG,yG])E_totalY = E_totalY*cut_total(charges,[xG,yG])ax.cla()# plot point chargesfor C in charges:    if C.q>0:        plot(C.pos[0], C.pos[1], 'bo', ms=8*sqrt(C.q))    if C.q<0:        plot(C.pos[0], C.pos[1], 'ro', ms=8*sqrt(-C.q))ax.quiver(xG,yG,E_totalX,E_totalY,pivot='middle',minshaft=0.1,minlength=0.3,headlength=2,headaxislength=2,headwidth=3,scale=4,alpha=0.4,width=0.002,linestyle='solid')xlabel('x')ylabel('y')ax.set_aspect(1)plt.savefig('fig2.png')#show()The output looks like follows:The problem with the first example is how to choose colormap. I played with different maps from http://matplotlib.org/examples/color/colormaps_reference.html but none gave me a really satisfying result. The problem with the first example is the cut function which doesn't work automatically. You need to cut out a specific region manually. So what would be your suggestions to improve those plots? If you have alternative FOSS tools, which can do those plots better or easier please share it. Please don't share LaTeX solutions since I asked a corresponding question on TeX.sx of how to attack this problem with PsTricks or something like that: https://tex.stackexchange.com/questions/225176/visualize-singular-vector-fields-with-tikz-or-pstricks-and-friendsI think the question should be community wiki. Edit Since especially the colormap part of the question is not answered yet, I add another simple example where the colormap doesn't look nice:%matplotlib inlinefrom pylab import *X=linspace(-2,2,40)Y=linspace(-2,2,24)X,Y=meshgrid(X, Y)def E(x,y):    r = sqrt(x**6 + y**6)    return (x/r,y/r)def E_dir(x,y):    #direction field    Ex,Ey=E(x,y)    n=sqrt(Ex**2+Ey**2)    return [Ex/n, Ey/n]Ex,Ey = E(X,Y)Exdir,Eydir = E_dir(X,Y)EE=sqrt(Ex**2+Ex**2)EQ  = quiver(X,Y,Exdir,Eydir,EE,cmap='autumn')show()Almost everything of the picture just looks red. 
Visually appealing ways to plot singular vector fields with matplotlib or other foss tools
python;data visualization
[I took your sample program as a starting point and adapted Colormap Normalization from the matplotlib wiki.]Almost everything of the picture just looks red.Indeed.  They problem is that there is a very narrow divergence in your data and because the colormap is scaled linearly almost all of the plot will be mapped to the lower limit of the colorbar.Q = quiver(X,Y,Exdir,Eydir,EE,cmap='autumn')That means we need to rescale the colorbar.  For a diveregence a logscale sounds good.  First of all we needimport matplotlib.colors as colorsso we can mess with the colors.  Next we simply putQ = quiver(X,Y,Exdir,Eydir,EE,cmap='autumn',           norm=colors.LogNorm(vmin=EE.min(),vmax=EE.max()))This will scale the colorbar logarithmically.That looks a little better already.  Now let's use another colormap because yellow is so hard to see.  I think coolwarm looks okay and suits the diverging data (I think Paraview uses this by default).Q = quiver(X,Y,Exdir,Eydir,EE,cmap='coolwarm',           norm=colors.LogNorm(vmin=EE.min(),vmax=EE.max()))As a last step let's add a nice colorbar indicating the divergence.fig = figure()Q = quiver(X,Y,Exdir,Eydir,EE,cmap='coolwarm',           norm=colors.LogNorm(vmin=EE.min(),vmax=EE.max()))fig.colorbar(Q,extend='max')Of course you can also put a colormap in the background as in the MATLAB answer.fig = figure()I = imshow(EE,extent=[np.min(X),np.max(X),np.min(Y),np.max(Y)],cmap='autumn')Q = quiver(X,Y,Exdir,Eydir)fig.colorbar(I,extend='max')or log scaledfig = figure()I = imshow(EE,extent=[np.min(X),np.max(X),np.min(Y),np.max(Y)],cmap='autumn',           norm=colors.LogNorm(vmin=EE.min(),vmax=EE.max()))Q = quiver(X,Y,Exdir,Eydir)fig.colorbar(I,extend='max')Another method is to clip the data to a maximum value.  The imshow command can do this on its own using the vmax key, for quiver you have to do it manually using numpy.clip.  This sets all values outside of a certain threshold to that threshold.numpy.clip(a, a_min, a_max, out=None)But because your function is changing so rapidly you probably have to use a logarithmic colormap in this case as well.fig = figure()Q = quiver(X,Y,Exdir,Eydir,clip(EE,None,1),cmap='autumn')fig.colorbar(Q,extend='max')