Inserting quality factors in rhone valley model
Hi @mafanasiev
I hope you are doing well! Would you have some time to help me over here with implementing 3D attenuation (via Federica's Lanza model) to the model for the Rhone valley? I think Federicas model is imposed on the same grid as Tobias Velocity model. Attached you can see the files (slightly different than the ones we used for reading Tobias model in the tools.py you implemented, but should have his model in the Vs, Vp).
I tried to write a function that reads these files and gets the Vs, Vp, Qs, Qp and insert it into the layered model:
- in visp_tools.py:
def read_veloq_model(filename: str) -> (xr.DataArray, xr.DataArray):
# """
# Read file that contains Q and V on a regular X-Y-Z grid, and
# return two DataArrays (V and Q) each with dims ("x", "y", "z").
# The file is assumed to have exactly this layout:
# • Line 1: some header text, e.g. "velocity X-Y-Z grid, Origin= … nzco= …"
# • Line 2: column names:
# latitude longitude X_km Y_km Z_km Q V Q-DWS Q-DRE
# • From line 3 onward: one row per grid point, with those nine columns in order.
# We:
# 1) np.loadtxt(skiprows=2) to grab all floating data.
# 2) Extract columns:
# – X_km (col 2), Y_km (col 3), Z_km (col 4)
# – Q (col 5), V (col 6)
# 3) Convert X_km→X_m, Y_km→Y_m, Z_km→Z_m (depth‐positive) exactly as read_velomod:
# X_m = X_km * 1000
# Y_m = Y_km * 1000
# Z_m = –Z_km * 1000 (so that Z_km=–5.00 → Z_m=+5000 m)
# 4) Find unique sorted coords in X_m, Y_m, Z_m and reshape Q/V into a 3D array
# shaped (n_z, n_y, n_x). Then transpose to (“x”,“y”,“z”) to match read_velomod.
# 5) Return two DataArrays: (1) velocity V(x,y,z) and (2) quality Q(x,y,z).
# Usage examples:
# vp_da, Qp_da = read_veloq_model("/…/vlxyzltln_Qpmodel.out")
# vs_da, Qs_da = read_veloq_model("/…/vlxyzltln_Qsmodel.out")
# # Or, equivalently (because velomod_P.out has a Q column):
# vp_da, Qp_da = read_veloq_model("/…/velomod_P.out")
# vs_da, Qs_da = read_veloq_model("/…/velomod_S.out")
# Returns
# -------
# velocity : xr.DataArray, dims=("x","y","z")
# quality factors : xr.DataArray, dims=("x","y","z")
# """
# 1) Load numeric data (skip the two‐line header)
data = np.loadtxt(filename, skiprows=2)
# Columns: [lat, lon, X_km(2), Y_km(3), Z_km(4), Q(5), V(6), Q-DWS(7), Q-DRE(8)]
# 2) Extract columns
X_km = data[:, 2]
Y_km = data[:, 3]
Z_km = data[:, 4]
Q_vals = data[:, 5]
V_vals = data[:, 6]
# 3) Convert to meters; note Z is inverted so that depth is positive
X_m = X_km * 1000.0
Y_m = Y_km * 1000.0
Z_m = -Z_km * 1000.0
# 4) Sort unique grid‐coordinates
x_coords = np.sort(np.unique(X_m))
y_coords = np.sort(np.unique(Y_m))
z_coords = np.sort(np.unique(Z_m))
nx, ny, nz = x_coords.size, y_coords.size, z_coords.size
# 5) Build empty 3D arrays for V and Q, shape (n_z, n_y, n_x)
V_3d = np.full((nz, ny, nx), np.nan, dtype=float)
Q_3d = np.full((nz, ny, nx), np.nan, dtype=float)
# 6) Make lookup maps from coordinate → index
xi_map = {x_coords[i]: i for i in range(nx)}
yi_map = {y_coords[i]: i for i in range(ny)}
zi_map = {z_coords[i]: i for i in range(nz)}
# 7) Fill the 3D arrays row by row
for row in data:
xm = row[2] * 1000.0
ym = row[3] * 1000.0
zm = -row[4] * 1000.0
qv = row[5]
vv = row[6]
i_x = xi_map[xm]
i_y = yi_map[ym]
i_z = zi_map[zm]
V_3d[i_z, i_y, i_x] = vv
Q_3d[i_z, i_y, i_x] = qv
# 8) Transpose to dims ("x","y","z") so that index order is (nx, ny, nz)
V_da = xr.DataArray(
data=np.transpose(V_3d, (2, 1, 0)), # (z,y,x)->(x,y,z)
coords={"x": x_coords, "y": y_coords, "z": z_coords},
dims=("x", "y", "z"),
)
Q_da = xr.DataArray(
data=np.transpose(Q_3d, (2, 1, 0)),
coords={"x": x_coords, "y": y_coords, "z": z_coords},
dims=("x", "y", "z"),
)
return V_da, Q_daI would like to ask you, does this function make sense? It seems to read the files and create the arrays for Vel and quality factors, but If you had some time to check it, whether it is similarly well parsing these values for the correct domain.
- Then in the snippet:
## Set up layers
surface = lm.interface.Surface(s_topo_gmrt)
bottom_of_roten = lm.interface.Surface(
s_topo_gmrt.assign_attrs(reference_elevation=lm.interface.Depth(560.0))
)
# vb = lm.material.elastic.Velocity.from_params(rho=vb_rho, vp=vb_vp, vs=vb_vs)
vb = lm.material.elastic.Velocity.from_params(rho=vb_rho, vp=vb_vp, vs=vb_vs)
bottom_of_vb = lm.interface.Hyperplane.at(-500.0)
# ─────────────────────────────────────────────────────────────────────────────
# First, read (vp, Qp) and (vs, Qs) from the Q‐model files and unpack them:
vp_da, Qp_da = visp_tools.read_veloq_model(
"/Users/maria/swiss_wavefields-main/data/vlxyzltln_Qpmodel.out"
)
vs_da, Qs_da = visp_tools.read_veloq_model(
"/Users/maria/swiss_wavefields-main/data/vlxyzltln_Qsmodel.out"
)
# ─────────────────────────────────────────────────────────────────────────────
# Now call Velocity.from_params exactly as before, passing only rho, vp, and vs:
bedrock = lm.material.elastic.Velocity.from_params(
rho = 2000,
vp = vp_da,
vs = vs_da,
qs = Qs_da,
qp = Qp_da,
)
print(bedrock)
bottom_of_mesh = lm.interface.Hyperplane.at(-45000.0)
# -
# ## Bring together in a layered model
layered_model = lm.LayeredModel(
[
surface,
lm.material.elastic.Velocity.from_dataset(
mod.isel(x=slice(None, None, 10), y=slice(None, None, 10))
),
bottom_of_roten,
vb,
bottom_of_vb,
bedrock,
bottom_of_mesh,
]
)
It doesn't work with the module:
lm.material.elastic.Velocity.from_params
giving: TypeError: from_params() got an unexpected keyword argument 'qs'
So a second question: What is the right module to read velocities, rho and quality factors? I checked the other modules (with_attenuation, visco.elastic.) but unsure how to properly use them.
Thanks a lot in advance for your time.
Best, maria
p.s.: the files where velocities and Qs are read from are attached now.