From 664a7bb484fd1823a651e48a361a40ec15ec53b2 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Mar 2016 13:26:50 -0800 Subject: [PATCH] Updated plotting functions for MT --- SimPEG/MT/Utils/dataUtils.py | 55 +++++++++++++++++++ SimPEG/Mesh/MeshIO.py | 45 ++++++++++++++- ...test_Problem1D_againstAnalyticHalfspace.py | 15 ----- 3 files changed, 99 insertions(+), 16 deletions(-) diff --git a/SimPEG/MT/Utils/dataUtils.py b/SimPEG/MT/Utils/dataUtils.py index b1947cb8..8df9d5d0 100644 --- a/SimPEG/MT/Utils/dataUtils.py +++ b/SimPEG/MT/Utils/dataUtils.py @@ -155,6 +155,61 @@ def plotMT1DModelData(problem,models,symList=None): ax.yaxis.set_tick_params(labelsize=fontSize) return fig +def plotImpAppRes(dataArrays,plotLoc,textStr=[]): + ''' Plots amplitude impedance and phase''' + # fig = plt.figure(1,(7, 7)) + import plotDataTypes as pDt + # axes = ImageGrid(fig, (0.05,0.05,0.875,0.875),nrows_ncols = (2, 2),axes_pad = 0.25,add_all=True,share_all=True,label_mode = "L") + # Make the figure and axes + fig,axT=plt.subplots(2,2,sharex=True) + axes = axT.ravel() + fig.set_size_inches((13.5,7.0)) + fig.suptitle('{:s}\nStation at: {:.1f}x ; {:.1f}y'.format(textStr,plotLoc[0],plotLoc[1])) + # Have to deal with axes + # Set log + for ax in axes.ravel(): + ax.set_xscale('log') + + axes[0].invert_xaxis() + axes[0].set_yscale('log') + axes[2].set_yscale('log') + # Set labels + axes[2].set_xlabel('Frequency [Hz]') + axes[3].set_xlabel('Frequency [Hz]') + axes[0].set_ylabel('Apperent resistivity [Ohm m]') + axes[1].set_ylabel('Apperent phase [degrees]') + axes[1].set_ylim(-180,180) + axes[2].set_ylabel('Impedance amplitude [V/A]') + axes[3].set_ylim(-180,180) + axes[3].set_ylabel('Impedance angle [degrees]') + + + # Plot the data + for nr,dataArray in enumerate(dataArrays): + if nr==1: + parSym = '*' + else: + parSym = 's' + # app res + pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zxy',par='res',pSym=parSym) + pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zyx',par='res',pSym=parSym) + # app phs + pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zxy',par='phs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zyx',par='phs',pSym=parSym) + # imp abs + pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxx',par='abs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxy',par='abs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyx',par='abs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyy',par='abs',pSym=parSym) + # imp abs + pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxx',par='phs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxy',par='phs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyx',par='phs',pSym=parSym) + pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyy',par='phs',pSym=parSym) + + return fig,axes + + def printTime(): import time print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()) diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py index 7501a66f..3ae4d88f 100644 --- a/SimPEG/Mesh/MeshIO.py +++ b/SimPEG/Mesh/MeshIO.py @@ -141,7 +141,6 @@ class TensorMeshIO(object): vtkObj.GetCellData().AddArray(vtkDoubleArr) # Set the active scalar vtkObj.GetCellData().SetActiveScalars(models.keys()[0]) - # vtkObj.Update() # Check the extension of the fileName ext = os.path.splitext(fileName)[1] @@ -158,6 +157,50 @@ class TensorMeshIO(object): vtrWriteFilter.SetFileName(fileName) vtrWriteFilter.Update() + def _toVTRObj(mesh,models=None): + """ + Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. + + Input: + :param str, path to the output vtk file + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells + + """ + # Import + from vtk import vtkRectilinearGrid as rectGrid, VTK_VERSION + from vtk.util.numpy_support import numpy_to_vtk + + # Deal with dimensionalities + if mesh.dim >= 1: + vX = mesh.vectorNx + xD = mesh.nNx + yD,zD = 1,1 + vY, vZ = np.array([0,0]) + if mesh.dim >= 2: + vY = mesh.vectorNy + yD = mesh.nNy + if mesh.dim == 3: + vZ = mesh.vectorNz + zD = mesh.nNz + # Use rectilinear VTK grid. + # Assign the spatial information. + vtkObj = rectGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + if models is not None: + for item in models.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + # Set the active scalar + vtkObj.GetCellData().SetActiveScalars(models.keys()[0]) + return vtkObj def readModelUBC(mesh, fileName): """ diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py index 66f5f41d..01373395 100644 --- a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py +++ b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py @@ -139,21 +139,6 @@ class TestAnalytics(unittest.TestCase): # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) - # def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr) - # def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp) - - # def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr) - # def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp) - - # def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr) - # def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp) - - # def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr) - # def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp) - - # def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr) - # def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp) - # Primary/secondary def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr) def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)