In [1]:

    

from os.path import join
import numpy as np

from scipy.interpolate import interp2d
from pyproj import Proj, transform

from mtpy.modeling.modem import Model, Data
from mtpy.utils import gis_tools

from mtpy.contrib.netcdf import nc

from mtpy.contrib.netcdf.modem_to_netCDF import interpolate, median_spacing


    

    




No handlers could be found for logger "mtpy.utils.decorator"
ERROR:root:the config yaml file /home/jami/.anaconda3/envs/mtpy_fork/lib/python2.7/site-packages/mtpy-1.0.1-py2.7.egg/mtpy/utils/logging.yml does not exist?
None
WARNING:mtpy.utils.decorator:GDAL_DATA environment variable is not set  Please see https://trac.osgeo.org/gdal/wiki/FAQInstallationAndBuilding#HowtosetGDAL_DATAvariable 
WARNING:mtpy.utils.decorator:GDAL_DATA environment variable is not set  Please see https://trac.osgeo.org/gdal/wiki/FAQInstallationAndBuilding#HowtosetGDAL_DATAvariable 






    




GDAL is not working !!!
GDAL is not working !!!
GDAL is not working !!!
If you want to write a vtk file for 3d viewing, you need download and install evtk from https://bitbucket.org/pauloh/pyevtk
Note: if you are using Windows you should build evtk first witheither MinGW or cygwin using the command: 
    python setup.py build -compiler=mingw32  or 
    python setup.py build -compiler=cygwin
If you want to write a vtk file for 3d viewing, you need to install pyevtk
Note: if you are using Windows you should build evtk first witheither MinGW or cygwin using the command: 
    python setup.py build -compiler=mingw32  or 
    python setup.py build -compiler=cygwin
If you want to write a vtk file for 3d viewing, you need download and install evtk from https://bitbucket.org/pauloh/pyevtk
If you want to write a vtk file for 3d viewing, you need to pip install PyEVTK: https://bitbucket.org/pauloh/pyevtk
Note: if you are using Windows you should build evtk first witheither MinGW or cygwin using the command: 
    python setup.py build -compiler=mingw32  or 
    python setup.py build -compiler=cygwin






    




If you want to write a vtk file for 3d viewing, you need download and install evtk from https://bitbucket.org/pauloh/pyevtk

In [2]:

    

# Define Data and Model Paths
MT_PATH = 'E:/Githubz/ummazannat/mtpy/examples/model_files/ModEM/'
# MT_PATH = "../../examples/model_files/ModEM/"

data = Data()
data.read_data_file(data_fn=join(MT_PATH, 'ModEM_Data.dat'))

# create a model object using the data object and read in model data
model = Model(data_obj=data)
model.read_model_file(model_fn=join(MT_PATH, 'ModEM_Model_File.rho'))


    

    




INFO:Data:Set rotation angle to 0.0 deg clockwise from N
/home/jami/.anaconda3/envs/mtpy_fork/lib/python2.7/site-packages/mtpy-1.0.1-py2.7.egg/mtpy/utils/calculator.py:322: RuntimeWarning: invalid value encountered in double_scalars
  z_rel_err = error/z_amp






    




      -15500.000      -21500.000         -67.000

    0.000

In [3]:

    

center = data.center_point

# dictionary for resistivity model data
# populate variables

resistivity_data = {
        'x': center.east.item() + (model.grid_east[1:] + model.grid_east[:-1])/
2,
        'y': center.north.item() + (model.grid_north[1:] + model.grid_north[:-1
])/2,
        'z': (model.grid_z[1:] + model.grid_z[:-1])/2,
        'resistivity': np.transpose(model.res_model, axes=(2, 0, 1))
    }

zone_number, is_northern, utm_zone = gis_tools.get_utm_zone(center.lat.item(), center.lon.item())
source_proj = Proj('+proj=utm +zone=%d +%s +datum=%s' % (zone_number, 'north' if is_northern else 'south', 'WGS84'))


    

In [4]:

    

# grid_proj is the lat lon global projection.
# At this moment it is hard coded because we are not using GDAL

grid_proj = Proj(init='epsg:4326')

result = interpolate(resistivity_data, source_proj, grid_proj,
                     center, median_spacing(model.grid_east), median_spacing(model.grid_north))

output_file = 'output.nc'

nc.write_resistivity_grid(output_file, grid_proj,
                          result['latitude'], result['longitude'], result['depth'],
                          result['resistivity'], z_label='depth')


    

In [ ]: