Skip to content

Commit

Permalink
Gmt Export (#238)
Browse files Browse the repository at this point in the history 
* Adding option to export output in text format

* Adding functionality to export MT resistivity from SGrid volumes
  • Loading branch information
@geojunky
geojunky committed Feb 18, 2022
1 parent b3fcde3 commit 7ebf1a7
Show file tree
Hide file tree
Showing 4 changed files with 276 additions and 48 deletions.
142 changes: 104 additions & 38 deletions seismic/receiver_fn/plot_ccp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
from PIL.PngImagePlugin import PngImageFile, PngInfo

from seismic.receiver_fn.rf_ccp_util import CCPVolume, CCP_VerticalProfile, CCP_DepthProfile, Gravity
from utils.sgrid import SGrid

CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
@click.group(context_settings=CONTEXT_SETTINGS)
Expand Down Expand Up @@ -93,7 +94,14 @@ def validate_start_end(line):
@click.argument('profile-def', type=click.Path(exists=True, dir_okay=False), required=True)
@click.option('--gravity-grid', type=click.Path(exists=True, dir_okay=False), default=None, show_default=True,
help='Gravity grid in tif/ers format. If provided, a gravity line-plot is added to each'
' vertical profile')
' vertical profile. Note that this parameter is ignored if the output-format is set '
'to txt')
@click.option('--mt-sgrid', type=click.Path(exists=True, dir_okay=False), default=None, show_default=True,
help='MT grid in sgrid format. If provided, an additional column with resistivity values are '
'added in the output text file. Note that this parameter is ignored if the output-format '
'is set to png')
@click.option('--mt-utm-zone', type=str, default=None, show_default=True,
help='UTM zone for mt-grid e.g. 53S. This parameter is ignored if mt-grid is not provided')
@click.option('--dx', type=float, default=5, show_default=True,
help='Horizontal distance-step (km) between start and end locations of a profile')
@click.option('--dz', type=float, default=0.5, show_default=True,
Expand Down Expand Up @@ -130,10 +138,11 @@ def validate_start_end(line):
help='Stations closer than max-station-dist (km) to a profile are labelled on plots')
@click.option('--output-folder', type=str, default='./', show_default=True,
help='Output folder')
@click.option('--output-format', type=str, default='png', show_default=True,
@click.option('--output-format', type=click.Choice(['png', 'txt']), default='png', show_default=True,
help='Output format')
def vertical(ccp_h5_volume, profile_def, gravity_grid, dx, dz, max_depth, swath_width, ds, extend, cell_radius,
idw_exponent, pw_exponent, amplitude_min, amplitude_max, max_station_dist, output_folder, output_format):
def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx, dz, max_depth, swath_width, ds,
extend, cell_radius, idw_exponent, pw_exponent, amplitude_min, amplitude_max, max_station_dist,
output_folder, output_format):
""" Plot CCP vertical profile \n
CCP_H5_VOLUME: Path to CCP volume in H5 format (output of rf_3dmigrate.py)\n
PROFILE_DEF: text file containing start and end locations of each vertical profile as\n
Expand All @@ -144,6 +153,10 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, dx, dz, max_depth, swath_
"""
log.setLevel(logging.DEBUG)

# sanity checks
if (mt_sgrid and not mt_utm_zone): assert 0, 'UTM zone for {} not specified. ' \
'See help for --mt-utm-zone'.format(mt_sgrid)

log.info('Loading CCP volume {}..'.format(ccp_h5_volume))
vol = CCPVolume(ccp_h5_volume)

Expand All @@ -156,6 +169,12 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, dx, dz, max_depth, swath_
gravity = Gravity(gravity_grid)
# end if

mt = None
if(mt_sgrid):
log.info('Loading MT SGrid {}..'.format(mt_sgrid))
mt = SGrid(mt_sgrid, mt_utm_zone)
# end if

for s, e in profiles:
log.info('Processing {}..'.format('-'.join([s, e])))

Expand All @@ -164,36 +183,69 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, dx, dz, max_depth, swath_
cell_radius=cell_radius, idw_exponent=idw_exponent,
pw_exponent=pw_exponent, max_station_dist=max_station_dist)

fig, gax, ax = None, None, None
if(gravity):
fig, (gax, ax) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 5]}, sharex=True)
fname = os.path.join(output_folder, '{}-{}.{}'.format(s, e, output_format))

if(output_format.lower() == 'txt'):
ggx, ggd = np.meshgrid(vprof._gx, vprof._gd)

grid = np.zeros((vprof._grid.shape[0], 5))
grid[:, 0] = vprof._grid[:, 1]
grid[:, 1] = vprof._grid[:, 2]
grid[:, 2] = ggx.flatten()
grid[:, 3] = ggd.flatten()
grid[:, 4] = vprof._grid_vals.flatten()

if(mt):
grid = np.hstack((grid, np.atleast_2d(np.zeros(grid.shape[0])).T))

grid[:, 5] = mt.get_values(grid[:, 0], grid[:, 1], grid[:, 3]*1e3, p=3)

"""
fig, ax = plt.subplots(1)
fig.set_size_inches(25,10)
fig.set_dpi(300)
vals = np.reshape(np.array(grid[:, 5]), (vprof._gd.shape[0], vprof._gx.shape[0]))
ax.pcolor(np.log1p(vals), cmap='seismic_r')
plt.savefig(fname+'.png')
"""

np.savetxt(fname, grid, header='lon lat distance(km) depth(km) amplitude(arb. units) resistivity(ohm m)',
fmt='%10.10f')
else:
np.savetxt(fname, grid, header='lon lat distance(km) depth(km) amplitude(arb. units)', fmt='%10.10f')
# end if
else:
fig, ax = plt.subplots()
# end if
fig, gax, ax = None, None, None
if(gravity):
fig, (gax, ax) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 5]}, sharex=True)
else:
fig, ax = plt.subplots()
# end if

fig.set_size_inches((20, 10))
fig.set_dpi(300)
fig.set_size_inches((20, 10))
fig.set_dpi(300)

# plot profile
vprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max, gax=gax, gravity=gravity)
# plot profile
vprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max, gax=gax, gravity=gravity)

fname = os.path.join(output_folder, '{}-{}.{}'.format(s, e, output_format))
fig.savefig(fname)

if(output_format.lower() == 'png'):
# write meta-data for translating digitization
d = {'lon1': vprof._lon1, 'lat1': vprof._lat1,
'lon2': vprof._lon2, 'lat2': vprof._lat2,
'x1':np.min(vprof._gx), 'y1':np.min(vprof._gd),
'x2':np.max(vprof._gx), 'y2':np.max(vprof._gd)}
sd = json.dumps(d)

img = PngImageFile(fname)
meta = PngInfo()
meta.add_text('profile_meta', sd)

img.save(fname, pnginfo=meta)
img.close()
fig.savefig(fname)

if(output_format.lower() == 'png'):
# write meta-data for translating digitization
d = {'lon1': vprof._lon1, 'lat1': vprof._lat1,
'lon2': vprof._lon2, 'lat2': vprof._lat2,
'x1':np.min(vprof._gx), 'y1':np.min(vprof._gd),
'x2':np.max(vprof._gx), 'y2':np.max(vprof._gd)}
sd = json.dumps(d)

img = PngImageFile(fname)
meta = PngInfo()
meta.add_text('profile_meta', sd)

img.save(fname, pnginfo=meta)
img.close()
# end if
# end if
# end for
# end
Expand Down Expand Up @@ -228,7 +280,7 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, dx, dz, max_depth, swath_
help='Maximum amplitude for colorbar normalization')
@click.option('--output-folder', type=str, default='./', show_default=True,
help='Output folder')
@click.option('--output-format', type=str, default='png', show_default=True,
@click.option('--output-format', type=click.Choice(['png', 'txt']), default='png', show_default=True,
help='Output format')
def depth(ccp_h5_volume, profile_def, dx, dy, dz, extend, cell_radius,
idw_exponent, pw_exponent, amplitude_min, amplitude_max, output_folder, output_format):
Expand Down Expand Up @@ -258,15 +310,29 @@ def depth(ccp_h5_volume, profile_def, dx, dy, dz, extend, cell_radius,
cell_radius=cell_radius, idw_exponent=idw_exponent,
pw_exponent=pw_exponent)

fig, ax = plt.subplots()
fig.set_size_inches((10, 10))
fig.set_dpi(300)
fname = os.path.join(output_folder, '{}-{}-depth-{}.{}'.format(s, e, int(depth), output_format))

# plot profile
dprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max)
if(output_format.lower() == 'txt'):
ggx, ggy = np.meshgrid(dprof._gx, dprof._gy)

fname = os.path.join(output_folder, '{}-{}-depth-{}.{}'.format(s, e, int(depth), output_format))
fig.savefig(fname)
grid = np.zeros((dprof._grid.shape[0], 5))
grid[:, 0] = dprof._grid[:, 1]
grid[:, 1] = dprof._grid[:, 2]
grid[:, 2] = ggx.flatten()
grid[:, 3] = ggy.flatten()
grid[:, 4] = dprof._grid_vals.flatten()

np.savetxt(fname, grid, header='lon lat x(km) y(km) amplitude(arb. units)', fmt='%10.10f')
else:
fig, ax = plt.subplots()
fig.set_size_inches((10, 10))
fig.set_dpi(300)

# plot profile
dprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max)

fig.savefig(fname)
# end if
# end for
# end

Expand Down
Empty file added utils/__init__.py
View file
Empty file.
11 changes: 1 addition & 10 deletions utils/resample_volume.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
from netCDF4 import Dataset
import click
import os.path as path
from seismic.xcorqc.utils import rtp2xyz

# define utility functions
def rtp2xyz(r, theta, phi):
Expand All @@ -34,16 +35,6 @@ def rtp2xyz(r, theta, phi):
return xout
# end func

def xyz2rtp(x, y, z):
rout = np.zeros((x.shape[0], 3))
tmp1 = x * x + y * y
tmp2 = tmp1 + z * z
rout[:, 0] = np.sqrt(tmp2)
rout[:, 1] = np.arctan2(np.sqrt(tmp1), z)
rout[:, 2] = np.arctan2(y, x)
return rout
# end func

class Resample:
def __init__(self, fname, lonf, latf, zf, nn=1, p=4):
"""
Expand Down

0 comments on commit 7ebf1a7

Please sign in to comment.