From ed72fba0639781f9fce9bdb7aa167597f692c86a Mon Sep 17 00:00:00 2001 From: GudniRos Date: Thu, 2 Jul 2015 14:00:37 -0700 Subject: [PATCH] Added EDI files read support. Fixed all srcMT to take 2 inputs. --- .travis.yml | 3 +- ...r2.ipynb => MT1Dinversion_Scipy2015.ipynb} | 7 +- simpegMT/DataMT.py | 55 +++++- simpegMT/ProblemMT1D/Problems.py | 7 + simpegMT/SurveyMT.py | 5 +- ...test_Problem1D_againstAnalyticHalfspace.py | 2 +- .../test_Problem1D_totalDvsPSvsAnalytic.py | 2 +- simpegMT/Utils/__init__.py | 1 + simpegMT/Utils/dataUtils.py | 7 +- simpegMT/Utils/ediFilesUtils.py | 171 ++++++++++++++++++ 10 files changed, 243 insertions(+), 17 deletions(-) rename notebooks/{MT1D inversion test-nr2.ipynb => MT1Dinversion_Scipy2015.ipynb} (99%) create mode 100644 simpegMT/Utils/ediFilesUtils.py diff --git a/.travis.yml b/.travis.yml index 316f710c..d9c2446e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,7 +29,8 @@ install: - cd simpegem/ - python setup.py build_ext --inplace - cd ../ - - echo export PYTHONPATH=$PYTHONPATH:/home/travis/build/simpegem/simpegem >> .bashrc- source .bashrc + - echo export PYTHONPATH=$PYTHONPATH:/home/travis/build/simpegem/simpegem >> .bashrc + - source .bashrc - cd simpegmt # Run test diff --git a/notebooks/MT1D inversion test-nr2.ipynb b/notebooks/MT1Dinversion_Scipy2015.ipynb similarity index 99% rename from notebooks/MT1D inversion test-nr2.ipynb rename to notebooks/MT1Dinversion_Scipy2015.ipynb index c617d865..7909e3b5 100644 --- a/notebooks/MT1D inversion test-nr2.ipynb +++ b/notebooks/MT1Dinversion_Scipy2015.ipynb @@ -66,12 +66,7 @@ " rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))\n", "# Source list\n", "srcList =[]\n", - "tD = False\n", - "if tD:\n", - " for freq in freqs:\n", - " srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))\n", - "else:\n", - " for freq in freqs:\n", + "for freq in freqs:\n", " srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma_0))\n", "# Make the survey\n", "survey = simpegmt.SurveyMT.SurveyMT(srcList)\n", diff --git a/simpegMT/DataMT.py b/simpegMT/DataMT.py index 19dd77d5..914d5909 100644 --- a/simpegMT/DataMT.py +++ b/simpegMT/DataMT.py @@ -1,8 +1,11 @@ from SimPEG import Survey, Utils, Problem, np, sp, mkvc +from simpegMT.Utils import rec2ndarr +import simpegMT from scipy.constants import mu_0 import sys from numpy.lib import recfunctions as recFunc + ############ ### Data ### ############ @@ -18,6 +21,11 @@ class DataMT(Survey.Data): # Pass the variables to the "parent" method Survey.Data.__init__(self, survey, v) + # # Import data + # @classmethod + # def fromEDIFiles(): + # pass + def toRecArray(self,returnType='RealImag'): ''' Function that returns a numpy.recarray for a SimpegMT impedance data object. @@ -26,8 +34,6 @@ class DataMT(Survey.Data): ''' - def rec2ndarr(x,dt=float): - return x.view((dt, len(x.dtype.names))) # Define the record fields dtRI = [('freq',float),('x',float),('y',float),('z',float),('zxxr',float),('zxxi',float),('zxyr',float),('zxyi',float),('zyxr',float),('zyxi',float),('zyyr',float),('zyyi',float)] dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex)] @@ -61,14 +67,55 @@ class DataMT(Survey.Data): if 'RealImag' in returnType: outArr = outTemp - if 'Complex' in returnType: + elif 'Complex' in returnType: # Add the real and imaginary to a complex number - outArr = np.empty(outTemp.shape,dtype=dtCP) for comp in ['freq','x','y','z']: outArr[comp] = outTemp[comp].copy() for comp in ['zxx','zxy','zyx','zyy']: outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy() + else: + raise NotImplementedError('{:s} is not implemented, as to be RealImag or Complex.') # Return return outArr + + @classmethod + def fromRecArray(cls, recArray, srcType='primary'): + """ + Class method that reads in a numpy record array to MTdata object. + + Only imports the impedance data. + + """ + if srcType=='primary': + src = simpegMT.SurveyMT.srcMT_polxy_1Dprimary + elif srcType=='total': + simpegMT.SurveyMT.srcMT_polxy_1DhomotD + else: + raise NotImplementedError('{:s} is not a valid source type for MTdata') + + # Find all the frequencies in recArray + uniFreq = np.unique(recArray['freq']) + srcList = [] + dataList = [] + for freq in uniFreq: + # Initiate rxList + rxList = [] + # Find that data for freq + dFreq = recArray[recArray['freq'] == freq] + # Find the impedance rxTypes in the recArray. + rxTypes = [ comp for comp in recArray.dtype.names if len(comp)==4 and 'z' in comp and 'r' in comp or 'i' in comp] + for rxType in rxTypes: + # Find index of not nan values in rxType + notNaNind = ~np.isnan(dFreq[rxType]) + + locs = rec2ndarr(dFreq[['x','y','z']][notNaNind].copy()) + rxList.append(simpegMT.SurveyMT.RxMT(locs,rxType)) + dataList.append(dFreq[rxType][notNaNind].data) + srcList.append(src(rxList,freq)) + + # Make a survey + survey = simpegMT.SurveyMT.SurveyMT(srcList) + dataVec = np.hstack(dataList) + return cls(survey,dataVec) \ No newline at end of file diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index ead5764d..156e6bb5 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -20,11 +20,18 @@ class eForm_psField(BaseMTProblem): # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. _fieldType = 'e_1d' _eqLocs = 'EF' + _sigmaPrimary = None def __init__(self, mesh, **kwargs): BaseMTProblem.__init__(self, mesh, **kwargs) self.fieldsPair = FieldsMT_1D + # self._sigmaPrimary = sigmaPrimary + + @property + def sigmaPrimary(self): + return self._sigmaPrimary + def getA(self, freq): """ diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 5f9b2c23..974a8a76 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -219,15 +219,16 @@ class srcMT_polxy_1Dprimary(srcMT): MT source for both polarizations (x and y) given a 1D primary models. It assigns fields calculated from the 1D model as fields in the full space of the problem. """ - def __init__(self, rxList, freq, sigma1d): + def __init__(self, rxList, freq): # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).' - self.sigma1d = sigma1d + self.sigma1d = None srcMT.__init__(self, rxList, freq) def ePrimary(self,problem): # Get primary fields for both polarizations + self.sigma1d = problem._sigmaPrimary eBG_bp = homo1DModelSource(problem.mesh,self.freq,self.sigma1d) return eBG_bp diff --git a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py index afb24d5d..79e3eb4a 100644 --- a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py +++ b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py @@ -63,7 +63,7 @@ def appRes_TotalFieldNorm(sigmaHalf): # Make the survey survey, sigma, mesh = setupSurvey(sigmaHalf) - problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh) + problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh,sigma) problem.pair(survey) # Get the fields diff --git a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py index 74af340e..428ee2c7 100644 --- a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py +++ b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py @@ -109,7 +109,7 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf): # Make the survey # Primary secondary surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,False) - problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh) + problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh,sigma) problemPS.pair(surveyPS) # Analytic data dataAna = calculateAnalyticSolution(surveyPS.srcList,mesh,sigma) diff --git a/simpegMT/Utils/__init__.py b/simpegMT/Utils/__init__.py index 0b98a3c3..b683f8b4 100644 --- a/simpegMT/Utils/__init__.py +++ b/simpegMT/Utils/__init__.py @@ -1,3 +1,4 @@ from MT1Dsolutions import * # Add the names of the functions from MT1Danalytic import * from dataUtils import * +from ediFilesUtils import * diff --git a/simpegMT/Utils/dataUtils.py b/simpegMT/Utils/dataUtils.py index 709de033..12d0090b 100644 --- a/simpegMT/Utils/dataUtils.py +++ b/simpegMT/Utils/dataUtils.py @@ -17,5 +17,8 @@ def getAppRes(MTdata): def appResPhs(freq,z): app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2 - app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi) - return app_res, app_phs \ No newline at end of file + app_phs = np.arctan2(z.imag,z.real)*(180/np.pi) + return app_res, app_phs + +def rec2ndarr(x,dt=float): + return x.view((dt, len(x.dtype.names))) \ No newline at end of file diff --git a/simpegMT/Utils/ediFilesUtils.py b/simpegMT/Utils/ediFilesUtils.py new file mode 100644 index 00000000..983ec192 --- /dev/null +++ b/simpegMT/Utils/ediFilesUtils.py @@ -0,0 +1,171 @@ +# Functions to import and export MT EDI files. +from SimPEG import mkvc +from scipy.constants import mu_0 +from numpy.lib import recfunctions as recFunc +import simpegMT +from simpegMT.Utils.dataUtils import rec2ndarr + +# Import modules +import numpy as np +import os, osr, sys, re + +class EDIimporter: + """ + A class to import EDIfiles. + + """ + _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit + _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI) + + # Properties + filesList = None + comps = None + + # Hidden properties + _outEPSG = None + _2out = None + + + def __init__(self, EDIfilesList, compList=None, outEPSG=None): + + # Set the fileList + self.filesList = EDIfilesList + # Set the components to import + if compList is None: + self.comps = ['ZXXR','ZXYR','ZYXR','ZYYR','ZXXI','ZXYI','ZYXI','ZYYI','ZXX.VAR','ZXY.VAR','ZYX.VAR','ZYY.VAR'] + else: + self.comps = compList + if outEPSG is not None: + self._outEPSG = outEPSG + + def __call__(self,comps=None): + + if comps is None: + return self._data + + return self._data[comps] + + def importFiles(self): + """ + Function to import EDI files into a object. + + + """ + + # Constants that are needed for convertion of units + + # Temp lists + tmpStaList = [] + + tmpCompList = ['freq','x','y','z'] + tmpCompList.extend(self.comps) + # Make the outarray + dtRI = [(compS.lower().replace('.',''),float) for compS in tmpCompList] + # Loop through all the files + for nrEDI, EDIfile in enumerate(self.filesList): + # Read the file into a list of the lines + with open(EDIfile,'r') as fid: + EDIlines = fid.readlines() + # Find the location + latD, longD, elevM = _findLatLong(EDIlines) + # Transfrom coordinates + transCoord = self._transfromPoints(longD,latD) + # Extract the name of the file (station) + EDIname = EDIfile.split(os.sep)[-1].split('.')[0] + # Arrange the data + staList = [EDIname, EDIfile, transCoord[0], transCoord[1], elevM[0]] + # Add to the station list + tmpStaList.extend(staList) + + # Read the frequency data + freq = _findEDIcomp('>FREQ',EDIlines) + # Make the temporary rec array. + tArrRec = ( np.nan*np.ones( (len(freq),len(dtRI)) ) ).view(dtRI) #np.concatenate((freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI) + # Add data to the array + tArrRec['freq'] = mkvc(freq,2) + tArrRec['x'] = mkvc(np.ones((len(freq),1))*transCoord[0],2) + tArrRec['y'] = mkvc(np.ones((len(freq),1))*transCoord[1],2) + tArrRec['z'] = mkvc(np.ones((len(freq),1))*elevM[0],2) + for comp in self.comps: + # Deal with converting units of the impedance tensor + if 'Z' in comp: + unitConvert = self._impUnitEDI2SI + else: + unitConvert = 1 + # Rotate the data since EDI x is *north, y *east but Simpeg uses x *east, y *north (* means internal reference frame) + key = [comp.lower().replace('.','').replace(s,t) for s,t in [['xx','yy'],['xy','yx'],['yx','xy'],['yy','xx']] if s in comp.lower()][0] + tArrRec[key] = mkvc(unitConvert*_findEDIcomp('>'+comp,EDIlines),2) + # Make a masked array + mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype) + try: + outTemp = recFunc.stack_arrays((outTemp,mArrRec)) + except NameError as e: + outTemp = mArrRec + + # Assign the data + self._data = outTemp + + # % Assign the data to the obj + # nOutData=length(obj.data); + # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data; + def _transfromPoints(self,longD,latD): + # Coordinates convertor + if self._2out is None: + src = osr.SpatialReference() + src.ImportFromEPSG(4326) + out = osr.SpatialReference() + if self._outEPSG is None: + # Find the UTM EPSG number + Nnr = 700 if latD < 0.0 else 600 + utmZ = int(1+(longD+180.0)/6.0) + self._outEPSG = 32000 + Nnr + utmZ + out.ImportFromEPSG(self._outEPSG) + self._2out = osr.CoordinateTransformation(src,out) + # Return the transfrom + return self._2out.TransformPoint(longD,latD) + +# Hidden functions +def _findLatLong(fileLines): + latDMS = np.array(fileLines[_findLine(' LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float) + longDMS = np.array(fileLines[_findLine(' LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float) + elevM = np.array([fileLines[_findLine(' ELEV=',fileLines)[0]].split('=')[1].split()[0]],float) + # Convert to D.ddddd values + latS = np.sign(latDMS[0]) + longS = np.sign(longDMS[0]) + latD = latDMS[0] + latS*latDMS[1]/60 + latS*latDMS[2]/3600 + longD = longDMS[0] + longS*longDMS[1]/60 + longS*longDMS[2]/3600 + return latD, longD, elevM + +def _findLine(comp,fileLines): + """ Find a line number in the file""" + # Line counter + c = 0 + # List of indices for found lines + found = [] + # Loop through all the lines + for line in fileLines: + if comp in line: + # Append if found + found.append(c) + # Increse the counter + c += 1 + # Return the found indices + return found + +def _findEDIcomp(comp,fileLines,dt=float): + """ + Extract the data vector. + + Returns a list of the data. + """ + # Find the data + headLine, indHead = [(st,nr) for nr,st in enumerate(fileLines) if re.search(comp,st)][0] + # Extract the data + nrVec = int(headLine.split()[-1]) + c = 0 + dataList = [] + while c < nrVec: + indHead += 1 + dataList.extend(fileLines[indHead].split()) + c = len(dataList) + return np.array(dataList,dt) \ No newline at end of file