0

I have file containing points under the columns "x-cord", "y-cord", "value". These are irregularly spaced. I am trying to make a contour plot of "value" and overlay this over the original domain. I gave up trying to do this in both pgfplots and matlab and thought I would give python a go. An answer in any of these scripts would be fine. The python script is as follows

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import numpy.ma as ma
from numpy.random import uniform, seed
from scipy.spatial import ConvexHull
#
# Loading data
filename = "strain.dat"
coordinates = []
x_c = []
y_c = []
z_c = []
xyz = open(filename)
title = xyz.readline()
for line in xyz:
    x,y,z = line.split()
    coordinates.append([float(x), float(y), float(z)])
    x_c.append([float(x)])
    y_c.append([float(y)])
    z_c.append([float(z)])
xyz.close()
#
# Rehaping and translating data
x_c=np.ravel(np.array(x_c))
y_c=np.ravel(np.array(y_c))
z_c=np.ravel(np.array(z_c))
x_c = x_c-100.0
y_c = y_c-100.0
#
# Checking the convex hull
points=np.column_stack((x_c,y_c))
hull = ConvexHull(points);
plt.plot(points[hull.vertices,0], points[hull.vertices,1], 'r--', lw=2)
plt.scatter(x_c, y_c, marker='o', s=5, zorder=10)
#
# Mapping the irregular data onto a regular grid and plotting
xic = np.linspace(min(x_c), max(x_c), 1000)
yic = np.linspace(min(y_c), max(y_c), 1000)
zic = griddata((x_c, y_c), z_c, (xic[None,:], yic[:,None]))
CS = plt.contour(xic,yic,zic,15,linewidths=0.5,colors='k')
CS = plt.contourf(xic,yic,zic,15,cmap=plt.cm.summer)
plt.colorbar() # draw colorbar
#
#plt.scatter(x_c, y_c, marker='o', s=5, zorder=10)
plt.axis('equal')
plt.savefig('foo.pdf', bbox_inches='tight')
plt.show()

and the output looks like

Plot of the contour plot overlaid with the original irregular points, the contour plot exceeds the edges of the original points due to the size of the convex hull

The problem is that griddata uses a convex hull and this convex hull exceeds the edges of the irregular data. Is there any way to set the values of the griddata points which are outside the edges of the boundary of the original points to zero?

Edit

In the end I threw in the towel and reverted back to Matlab. I'll have to export the data to pgfplots to get a nice plot. The code I came up with was 

x = strain.x;
y = strain.y;
z = strain.eps;
% Get the alpha shape (couldn't do this in python easily)
shp = alphaShape(x,y,.001);
% Get the boundary nodes
[bi, xy] = boundaryFacets(shp);
no_grid = 500;
xb=xy(:,1);
yb=xy(:,2);
[X,Y] = ndgrid(linspace(min(x),max(x),no_grid),linspace(min(y),max(y),no_grid));
Z = griddata(x,y,z,X,Y,'v4');

% Got through the regular grid and set the values which are outside the boundary of the original domain to Nans
for j = 1:no_grid
    [in,on] = inpolygon(X(:,j),Y(:,j),xb,yb);
    Z(~in,j) = NaN;
end

contourf(X,Y,Z,10),axis equal
colorbar
hold on
plot(xb,yb)
axis equal
hold off

Here is the resulting image.enter image description here

If someone can do something similar in Python I'll happily accept the answer.

1QuickQuestion
  • 417
  • 6
  • 16
  • 2
    For a nice trapezoid, each *row* has two points - a min and a max; you could compare the points in the griddata to the boundary *by row* and determine if they are between (True) or outside (False) -creating a boolean mask then multiply the griddata by the mask. – wwii Jul 11 '17 at 14:17
  • That is a nice idea, it's a bit beyond my python capabilities though as the irregular points are not on nice regular rows. – 1QuickQuestion Jul 11 '17 at 15:37
  • In order for someone to give you a solution, you need to provide a [mcve] of the problem. @wwii's comment is in principle correct, so it all depends on how the data looks like how easy or not this is. – ImportanceOfBeingErnest Jul 11 '17 at 15:54

1 Answers1

0

I had to plot interpolated data on a complex geometry (see the blue points on figure) P(x,z) (z is the horizontal coordinate). I used mask operations and it worked well. Without mask, the whole square (x=0..1 ; z=0..17.28) is covered by contourf.

## limiting values for geometry
xmax1=0.408
zmin1=6.
xmax2=0.064
zmin2=13.12
xmin=0.
xmax=1.
zmin=0.
zmax=17.28

# Grid for points
x1 = np.arange(xmin,xmax+dx,dx)
z1 = np.arange(zmin,zmax+dz,dz)
zi2,xi2 = np.meshgrid(z1,x1)
mask = (((zi2 > zmin2) & (xi2 > xmax2)) | ((zi2 > zmin1) & (zi2 <= zmin2) & (xi2 > xmax1)))
zim=np.ma.masked_array(zi2,mask)
xim=np.ma.masked_array(xi2,mask)

# Grid for P values
# npz=z coordinates of data, npx is the x coordinates and npp is P values
grid_p = scipy.interpolate.griddata((npz, npx), npp, (zim,xim),method='nearest')
pm=np.ma.masked_array(grid_p,mask)

# plot 
plt.contour(zim, xim, pm, 25, linewidths=0.5, colors='k',corner_mask=False)
plt.contourf(zim, xim, pm, 25,vmax=grid_p.max(), vmin=grid_p.min(),corner_mask=False)
plt.colorbar()

# Scatter plot to check 
plt.scatter(npz,npr, marker='x', s=2)

plt.show()

enter image description here

alain
  • 11
  • 1