Depending on the parameter passing through the parameter NetCDF file, output data can be embedded in TIFF files or 5D NetCDF containers. We focus here on the treatment of NetCDF output data.
Data could be quickly viewed using ncview, paraview, or imageJ with the dedicated plugin.
For data post-treatment python could be used:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from Scientific.N import *
from Scientific.IO.NetCDF import *
from pylab import *
font = {'family':'serif','serif':'serif','weight':'normal','size':'14'}
mpl.rc('font',**font)
mpl.rc('text',usetex=True)
#-------------------------------------------------------
#Import data
#-------------------------------------------------------
def loadNC(path,var_ind,time0=0,time1=0):
file = NetCDFFile(path,"r")
#i = variable number / variables are stored in inverse order (e.g 'var_0 = W', 'var_1 = V', 'var_2 = U', with U,V,W the displacement fields along x,y and Z directions)
dimNames = file.dimensions.keys()
dimVals = file.dimensions.values()
var_name = file.variables.keys()
tmp = file.variables[var_name[var_ind]]
varShape = tmp.shape
# By default the 5 dimensions exist and some of them can be empty. Dimensions are stored as following : [t,c,z,y,x]
# - x,y, z are spatial dimensions (z is empty if 2D case, and y is empty if 1D case)
# - c (color channel) is never used (empty)
# - t (time) corresponds to the sequence of treated image (empty if only one couple of picture)
# To remove empty dimensions the following
var = np.squeeze(tmp[time0:time1,:,:,:,:])
file.close()
return var
#-------------------------------------------------------
#Strain calcul
#-------------------------------------------------------
#win = Correlation ZOI size (assumed squared)
#U,V or W = Var[i]
def strain(gridX,gridY,U,V,gridZ=0,W=np.zeros([1])):
E = {}
if len(W) == 0:
E['xx'] = np.zeros([U.shape[0],U.shape[1]])
E['yy'] = np.zeros([U.shape[0],U.shape[1]])
E['xy'] = np.zeros([U.shape[0],U.shape[1]])
E['yx'] = np.zeros([U.shape[0],U.shape[1]])
[E['xy'],E['xx']] = np.gradient(U,gridX,gridY)
[E['yy'],E['yx']] = np.gradient(V,gridX,gridY)
else:
E['xx'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['yy'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['zz'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['xy'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['xz'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['yx'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['yz'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['zx'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
E['zy'] = np.zeros([U.shape[0],U.shape[1],U.shape[2]])
[E['xy'],E['xx'],E['xz']] = np.gradient(U,gridX,gridY,gridZ)
[E['yy'],E['yx'],E['yz']] = np.gradient(V,gridX,gridY,gridZ)
[E['zy'],E['zx'],E['zz']] = np.gradient(W,gridX,gridY,gridZ)
return E
#-------------------------------------------------------
#Plot 2D
#-------------------------------------------------------
# Displacement field ("mode.nc") superimposed on picture ("deform.nc"), using grid node ("node.nc")
def plot2d(picture, field, gridX, gridY, title):
mini = field[:,:].mean() - 3*field[:,:].std()
maxi = field[:,:].mean() + 3*field[:,:].std()
fig = plt.figure()
plt.imshow(picture, origin='upper', interpolation='bicubic', cmap=cm.gray)
ax = plt.imshow(field, origin='upper', interpolation='bicubic', extent=(gridX.min(), gridX.max(), gridY.min(), gridY.max()), cmap=cm.jet, alpha=0.75)
fig.colorbar(ax, orientation='vertical')
plt.clim([mini, maxi])
plt.title(title)
plt.xticks(())
plt.yticks(())
plt.show()
plt.clf()
#-------------------------------------------------------
#--------------------PROGRAM----------------------------
time = 1
#modes
path = './modes.nc'
W = loadNC(path,2,time,time+1)
V = loadNC(path,1,time,time+1)
U = loadNC(path,0,time,time+1)
#deformed deformed image
path = './deform.nc'
images = loadNC(path,0,0,1)
#or (warning: for 3D Tiff use import tifffile
#http://www.lfd.uci.edu/~gohlke/code/tifffile.py.html and use it #as follows: name = './image.tif'
#with tifffile.TIFFfile(name) as imfile:
# Img = double(imfile.asarray())
#)
image = plt.imread('frame00000.tiff')
#grid
path = './Correlation/Limodin11_OFFEM_nodes.nc'
gridZ = loadNC(path,2,0,1)
gridY = loadNC(path,1,0,1)
gridX = loadNC(path,0,0,1)
#gradient step
dx = np.gradient(gridX)[0]
dy = np.gradient(gridY)[1]
dz = np.gradient(gridZ)[2]
x = (gridX).astype(int16)
y = (gridY).astype(int16)
z = (gridZ).astype(int16)
path = 'test_mask_modes.nc'
msk = loadNC(path,0,0,1)
Umask = np.ma.masked_array(U,np.abs(msk-1))
Vmask = np.ma.masked_array(V,np.abs(msk-1))
Wmask = np.ma.masked_array(W,np.abs(msk-1))
#strain
E = strain(dx,dy,Umask,Vmask,dz,Wmask)
plot2d(image, E['yy'][0,:,:], gridX[0,:,:], gridY[0,:,:], 'axial strain')
#-------------------------------------------------------
#-------------------------------------------------------
