Commit eaf621d4 authored by Mijian Xu's avatar Mijian Xu 😷

add dt2xa, dt2ya

parent 671b53c9
Pipeline #740 passed with stage
in 1 minute and 47 seconds
...@@ -16,9 +16,12 @@ class Client(Cmd): ...@@ -16,9 +16,12 @@ class Client(Cmd):
self.prompt = 'MGPro>' self.prompt = 'MGPro>'
self.aliases = {'c': self.do_continuation, self.aliases = {'c': self.do_continuation,
'g': self.do_gradient, 'g': self.do_gradient,
'p': self.do_plot,
'q': self.do_quit, 'q': self.do_quit,
'r': self.do_read, 'r': self.do_read,
'w': self.do_write, 'w': self.do_write,
'x': self.do_dt2xa,
'y': self.do_dt2ya,
'z': self.do_dt2za, 'z': self.do_dt2za,
'h': self.do_help} 'h': self.do_help}
self.result = None self.result = None
...@@ -93,9 +96,46 @@ horizontal gradient in azimuth of 0, 45, 90 or 135 or module of horizontal gradi ...@@ -93,9 +96,46 @@ horizontal gradient in azimuth of 0, 45, 90 or 135 or module of horizontal gradi
self.result = self.mg.gradient(arg) self.result = self.mg.gradient(arg)
def do_dt2za(self, arg): def do_dt2za(self, arg):
if len(arg.split()) != 2: '''
dt2za (z): Vertical magnetization
Syntax: dt2za i0 d0
i0: Magnetic dip of the region
d0: Magnetic declination of the region
'''
if len(arg.split()) == 2:
i0, d0 = [float(value) for value in arg.split()] i0, d0 = [float(value) for value in arg.split()]
self.result = self.mg.dt2za(i0, d0) self.result = self.mg.dt2za(i0, d0)
else:
print('Two argumrnts required')
def do_dt2xa(self, arg):
'''
dt2xa (x): Horizontal magnetization along x direction
Syntax: dt2xa i0 d0
i0: Magnetic dip of the region
d0: Magnetic declination of the region
'''
if len(arg.split()) == 2:
i0, d0 = [float(value) for value in arg.split()]
self.result = self.mg.dt2xa(i0, d0)
else:
print('Two argumrnts required')
def do_dt2ya(self, arg):
'''
dt2xa (y): Horizontal magnetization along y direction
Syntax: dt2ya i0 d0
i0: Magnetic dip of the region
d0: Magnetic declination of the region
'''
if len(arg.split()) == 2:
i0, d0 = [float(value) for value in arg.split()]
self.result = self.mg.dt2ya(i0, d0)
else:
print('Two argumrnts required')
def do_write(self, arg): def do_write(self, arg):
''' '''
...@@ -125,6 +165,14 @@ Syntax: write filename [\'lalo\'] ...@@ -125,6 +165,14 @@ Syntax: write filename [\'lalo\']
print('{}'.format(e)) print('{}'.format(e))
return return
def do_plot(self, arg):
'''
plot (p): Preview the result.
Syntax: plot
'''
self.mg.pltmap(self.result)
def completedefault(self, *args): def completedefault(self, *args):
readline.set_completer_delims(' \t\n;') readline.set_completer_delims(' \t\n;')
readline.parse_and_bind("tab: complete") readline.parse_and_bind("tab: complete")
......
import numpy as np import numpy as np
from mgpro import expand # from mgpro import expand
from mgpro.proj import * import expand
# from mgpro.proj import *
from proj import *
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.colors as colors import matplotlib.colors as colors
from scipy.fftpack import fft2, fftshift, ifft2, ifftshift from scipy.fftpack import fft2, fftshift, ifft2, ifftshift
from scipy.interpolate import griddata from scipy.interpolate import griddata
import argparse import argparse
np.warnings.filterwarnings('ignore')
def msg(): def msg():
return ''' return '''
...@@ -42,6 +45,12 @@ def xyz2grd(raw_data, dx, dy, xy_limit=None): ...@@ -42,6 +45,12 @@ def xyz2grd(raw_data, dx, dy, xy_limit=None):
grid_data = griddata(data[:, 0:2], data[:, 2], (ymesh, xmesh)) grid_data = griddata(data[:, 0:2], data[:, 2], (ymesh, xmesh))
return xrange, yrange, grid_data return xrange, yrange, grid_data
def direction_cos(i0, d0):
P0 = np.cos(np.deg2rad(i0)) * np.cos(np.deg2rad(d0))
Q0 = np.cos(np.deg2rad(i0)) * np.sin(np.deg2rad(d0))
R0 = np.sin(np.deg2rad(i0))
return P0, Q0, R0
def norm_uv(data, dx, dy): def norm_uv(data, dx, dy):
(len_row, len_col) = data.shape (len_row, len_col) = data.shape
...@@ -125,15 +134,33 @@ class mgmat(object): ...@@ -125,15 +134,33 @@ class mgmat(object):
return grad return grad
def dt2za(self, i0, d0): def dt2za(self, i0, d0):
P0 = np.cos(np.deg2rad(i0)) * np.cos(np.deg2rad(d0)) P0, Q0, R0 = direction_cos(i0, d0)
Q0 = np.cos(np.deg2rad(i0)) * np.sin(np.deg2rad(d0))
R0 = np.sin(np.deg2rad(i0))
u, v = norm_uv(self.data_sf, self.dx, self.dy) u, v = norm_uv(self.data_sf, self.dx, self.dy)
phi = np.sqrt(u**2 + v**2) / ((P0*u + Q0*v)*1j + R0*np.sqrt(u**2 + v**2)) phi = np.sqrt(u**2 + v**2) / ((P0*u + Q0*v)*1j + R0*np.sqrt(u**2 + v**2))
phi[np.where(np.isnan(phi))] = 0 + 0j
result = np.real(ifft2(ifftshift(phi * self.data_sf)))[self.row_begin: self.row_end + 1, self.col_begin: self.col_end + 1] result = np.real(ifft2(ifftshift(phi * self.data_sf)))[self.row_begin: self.row_end + 1, self.col_begin: self.col_end + 1]
return result return result
def pltmap(self, fig, data, breakpoint=None): def dt2xa(self, i0, d0):
P0, Q0, R0 = direction_cos(i0, d0)
u, v = norm_uv(self.data_sf, self.dx, self.dy)
trans_func = u*1j / ((P0*u + Q0*v)*1j + R0*np.sqrt(u**2 + v**2))
trans_func[np.where(np.isnan(trans_func))] = 0 + 0j
result = np.real(ifft2(ifftshift(trans_func * self.data_sf)))[self.row_begin: self.row_end + 1, self.col_begin: self.col_end + 1]
return result
def dt2ya(self, i0, d0):
P0, Q0, R0 = direction_cos(i0, d0)
u, v = norm_uv(self.data_sf, self.dx, self.dy)
trans_func = v*1j / ((P0*u + Q0*v)*1j + R0*np.sqrt(u**2 + v**2))
trans_func[np.where(np.isnan(trans_func))] = 0 + 0j
result = np.real(ifft2(ifftshift(trans_func * self.data_sf)))[self.row_begin: self.row_end + 1, self.col_begin: self.col_end + 1]
return result
def pltmap(self, data, breakpoint=None):
fig = plt.figure()
f_width = fig.get_figwidth() f_width = fig.get_figwidth()
f_height = fig.get_figheight() f_height = fig.get_figheight()
frac_w = min([f_height, f_width])/max([f_height, f_width]) frac_w = min([f_height, f_width])/max([f_height, f_width])
...@@ -141,12 +168,13 @@ class mgmat(object): ...@@ -141,12 +168,13 @@ class mgmat(object):
real_width *= frac_w real_width *= frac_w
fig.clf() fig.clf()
ax_raw = fig.gca() ax_raw = fig.gca()
pcm = ax_raw.pcolor(data, cmap='jet') pcm = ax_raw.imshow(data, cmap='jet', extent=[self.x[0], self.x[-1], self.y[0], self.y[-1]])
# ax_raw.figure.canvas.draw() # ax_raw.figure.canvas.draw()
cb = fig.colorbar(pcm, extend='both') cb = fig.colorbar(pcm, extend='both')
# ax_raw.set_position([.1, .125, real_width, real_height], which='original') # ax_raw.set_position([.1, .125, real_width, real_height], which='original')
# cb.ax.set_position([.8, .1, real_height/6.27, real_height]) # cb.ax.set_position([.8, .1, real_height/6.27, real_height])
fig.canvas.draw() fig.canvas.draw()
plt.show()
def savetxt(self, filename, result, to_latlon=False): def savetxt(self, filename, result, to_latlon=False):
points = np.zeros([result.size, 3]) points = np.zeros([result.size, 3])
...@@ -164,6 +192,4 @@ class mgmat(object): ...@@ -164,6 +192,4 @@ class mgmat(object):
if __name__ == '__main__': if __name__ == '__main__':
mg = mgmat('/Users/xumj/Codes/MGPro/example/mag_proj.dat', 2000, 2000) mg = mgmat('/Users/xumj/Codes/MGPro/example/mag_proj.dat', 2000, 2000)
mg.dt2za(46, -1) mg.pltmap(mg.dt2ya(46, -1))
# exec()
...@@ -5,7 +5,7 @@ packages = find_packages() ...@@ -5,7 +5,7 @@ packages = find_packages()
with open("README.md", "r") as fh: with open("README.md", "r") as fh:
long_description = fh.read() long_description = fh.read()
VERSION = "0.1.4" VERSION = "0.1.5"
setup(name='mgpro', setup(name='mgpro',
version=VERSION, version=VERSION,
author='Mijian Xu', author='Mijian Xu',
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment