Added EDI files read support.

Fixed all srcMT to take 2 inputs.

.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

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",

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)

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):
         """

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
 

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

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)

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 *

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
+    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)))

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)