<jupyter_start><jupyter_text>Load GBIS inversion result into HDF5 format<jupyter_code>%matplotlib inline
import os
import numpy as np
import pyproj
import utm
from matplotlib import pyplot as plt
from mintpy.utils import readfile
from mintpy import load_gbis, view

model_dir = os.path.expanduser('~/Papers/2021_Kirishima/figs_src/model')
os.chdir(model_dir)
print('Go to directory', model_dir)<jupyter_output>Go to directory /Users/yunjunz/Papers/2021_Kirishima/figs_src/model<jupyter_text>load GBIS result into HDF5<jupyter_code>load_gbis.main(['Shinmoe2008post/invert_1_2_C/invert_1_2_C.mat'])
load_gbis.main(['Shinmoe2017pre/invert_1_2_C/invert_1_2_C.mat'])
load_gbis.main(['Iwo2017pre/invert_1_2_C/invert_1_2_C.mat'])
load_gbis.main(['Iwo2017post/invert_1_2_C_C2/invert_1_2_C_C2.mat'])<jupyter_output>read mat file: /Users/yunjunz/Papers/2021_Kirishima/figs_src/model/Shinmoe2008post/invert_1_2_C/invert_1_2_C.mat
number of output HDF5 file: 2
creating figure in size of [12, 6]
------------------------------
read mask from file: /Users/yunjunz/Papers/2021_Kirishima/figs_src/model/data/KirishimaAlosAT424_20080929_20100705.mat
delete exsited file: /Users/yunjunz/Papers/2021_Kirishima/figs_src/model/Shinmoe2008post/invert_1_2_C/KirishimaAlosAT424_20080929_20100705.h5
create HDF5 file: /Users/yunjunz/Papers/2021_Kirishima/figs_src/model/Shinmoe2008post/invert_1_2_C/KirishimaAlosAT424_20080929_20100705.h5 with w mode
create dataset /hgt      of float32    in size of (361, 361)           with compression=None
create dataset /data     of float32    in size of (361, 361)           with compression=None
create dataset /model    of float32    in size of (361, 361)           with compression=None
create dataset /residual of float32    in size of (361, 361)           with compression=None
finishe[...]<jupyter_text>Backup 1: Plot observations with height below / above the free surface - Iwo-yama<jupyter_code>fname1 = os.path.join(model_dir, 'Iwo2017pre/invert_1_2_C/KirishimaAlos2AT131_20150106_20171010.h5')
fname2 = os.path.join(model_dir, 'Iwo2017pre/invert_1_2_C/KirishimaAlos2DT23_20150209_20170918.h5')

# Figure 2 bounding box in X/Y
S, N, W, E = 31.939, 31.954, 130.843, 130.861
ref_lat, ref_lon = 31.947, 130.853
ref_x, ref_y = utm.from_latlon(ref_lat, ref_lon)[:2]
pts = [[N, W], [N, E], [S, E], [S, W]]
xs = []
ys = []
for pt in pts:
    x, y = utm.from_latlon(pt[0], pt[1])[:2]
    xs.append(x - ref_x)
    ys.append(y - ref_y)
print('X min/max: {:.1f} / {:.1f}'.format(min(xs), max(xs)))
print('Y min/max: {:.1f} / {:.1f}'.format(min(ys), max(ys)))

# view options
cmd = f'view.py {fname1} data --noreference --noverbose --figtitle Iwo-yama --nocbar --notick '
cmd += ' --sub-lat {} {} --sub-lon {} {} '.format(S, N, W, E)   # Fig. 2
#cmd += ' --sub-lat 31.935 31.96  --sub-lon 130.832 130.870 '  # GBIS bounding box

atr, inps = view.prep_slice(cmd)[1:]
extent = (inps.geo_box[0], inps.geo_box[2], inps.geo_box[3], inps.geo_box[1])
dem   = readfile.read(fname1, datasetName='hgt',  box=inps.pix_box)[0]

# plot
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=[12, 5])
for i, (ax, fname, label) in enumerate(zip(axs, [fname1, fname2], ['asc', 'desc'])):
    data = readfile.read(fname, datasetName='data', box=inps.pix_box)[0]
    view.plot_slice(ax, data, atr, inps)
    ax.imshow(dem < 1300, vmin=0, vmax=1, cmap='gray_r', origin='upper', extent=extent, alpha=0.2, zorder=2)
    ax.set_title(label, fontsize=16)
plt.show()<jupyter_output>X min/max: -958.5 / 771.5
Y min/max: -903.1 / 789.0<jupyter_text>Backup 2: Plot observations below / above the free surface - Shinmoe-dake<jupyter_code>fname1 = os.path.join(model_dir, 'Shinmoe2008post/invert_1_2_C/KirishimaAlosAT424_20080929_20100705.h5')
fname2 = os.path.join(model_dir, 'Shinmoe2008post/invert_1_2_C/KirishimaAlosDT73_20081012_20100302.h5')

# Figure 2 bounding box in X/Y
S, N, W, E = 31.8975, 31.925, 130.868, 130.902
ref_lat, ref_lon = 31.91, 130.88
ref_x, ref_y = utm.from_latlon(ref_lat, ref_lon)[:2]
pts = [[N, W], [N, E], [S, E], [S, W]]
xs = []
ys = []
for pt in pts:
    x, y = utm.from_latlon(pt[0], pt[1])[:2]
    xs.append(x - ref_x)
    ys.append(y - ref_y)
print('X min/max: {:.1f} / {:.1f}'.format(min(xs), max(xs)))
print('Y min/max: {:.1f} / {:.1f}'.format(min(ys), max(ys)))

# view options
cmd = f'view.py {fname1} data --noreference --noverbose --figtitle Shinmoe-dake --nocbar --notick '
cmd += ' --sub-lat 31.8975 31.925 --sub-lon 130.868 130.902 '  # Fig. 2
#cmd += ' --sub-lat 31.935 31.96  --sub-lon 130.832 130.870 '  # GBIS bounding box

atr, inps = view.prep_slice(cmd)[1:]
extent = (inps.geo_box[0], inps.geo_box[2], inps.geo_box[3], inps.geo_box[1])
dem = readfile.read(fname1, datasetName='hgt', box=inps.pix_box)[0]

# plot
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=[12, 5])
for i, (ax, fname, label) in enumerate(zip(axs, [fname1, fname2], ['asc', 'desc'])):
    data = readfile.read(fname, datasetName='data', box=inps.pix_box)[0]
    view.plot_slice(ax, data, atr, inps)
    ax.imshow(dem < 1100, vmin=0, vmax=1, cmap='gray_r', origin='upper', extent=extent, alpha=0.3, zorder=2)
    ax.set_title(label, fontsize=16)
plt.show()<jupyter_output>X min/max: -1163.5 / 2104.9
Y min/max: -1405.5 / 1699.3