UCVM shakeout2 velocity model: Difference between revisions

Latest revision as of 00:38, 6 October 2026

Model overview

The shakeout2 velocity model for the ShakeOut Scenario 2.0 is based on CVM-S4.26.M01(with embedded GTL layer for basin regions) and GTL Tapering for all non-basin area.

Tapering's z range is from 0m to 700m, and material property floors used are vs=200m, vp=500m and density=200m a customized z coordinates list is pre-defined

four corners:


Mesh dimensions: 516 x 468 x 152

Grid 4 corners:
  -120.028425 31.874269
  -114.584135 31.887394
  -120.181534 36.079919
  -114.461905 36.095287

Creating the Data files with ucvm2mesh

UCVM's ucvm2mesh is used to generate 3 binary mesh data files: vs, vp and density for ROI

In 'debug mode', 3 coordinate lists are created: one for the lats, the lons and the z depths.

With these data files and lists, a simple NetCDF file is created

ucvm2mesh configuration file:


# List of CVMs to query
ucvmlist=cvmsi,elygtl:taper

# (optional) zranges: -Z zmin,zmax 
min_zrange=0.0
max_zrange=700.0

# (optional) floors: -L vs_floor,vp_floor,density_floor
vs_floor=200.0
vp_floor=500.0
density_floor=200.0

# UCVM conf file
ucvmconf=YOUR_UCVM_INSTALL_PATH/conf/ucvm.conf

# Gridding cell: CENTER or VERTEX 
gridtype=CENTER

# Spacing of cells (m)
spacing=1000.0
# (optional) Spacing of z depth (m)
# if not set, z_spacing = spacing 
#z_spacing = 500.0

# (optional) Allow varying z depth
# format: 1 row 1 target depth 
z_file = so2_zlist

# (optional) output latitudes filename
y_file = so2_ylist

# (optional) output longitudes filename
x_file = so2_xlist

# Projection
proj=+proj=utm +datum=WGS84 +zone=11
#to rotate, rot=-39.9
rot=0

# Origin 
# x,y: latlon
# z: depth (m)
x0=-120.033556
y0=31.869638
z0=0.0

# Number of cells along each dim
nx=516
ny=468
nz=152
# how to determine distance in meter for each dim
# x = nx * spacing
# y = ny * spacing
# z = nz * (spacing or z_spacing), or = z_file[idx])

# Partitioning of grid among processors
# relevant to chunking, should be divisible into nx,ny,nz
px=2
py=2
pz=2

# Vs/Vp minimum
vp_min=0
vs_min=0

# Mesh grid files
meshfile=so2.media
gridfile=so2.grid

# Mesh Grid formats
# IJK-12, IJK-20, IJK-32 : single triplet file
# SORD : separate vs, vp, rho files
meshtype=SORD

# Location of scratch dir
scratch=./scratch

ucvm2mesh z depths list: so2_zlist

latitude coordinate list: so2_xlist

longitude coordinate list: so2_ylist

ucvm2mesh command used:

ucvm2mesh -f so2.conf > so2.out 2>&1

and if so2.grid is known,

ucvm2mesh -f so2.conf -g so2.grid > so2.out 2>&1

Create NetCDF mesh data

make_nc.py

  
#!/usr/bin/env python3

import numpy as np
from netCDF4 import Dataset

## Update filenames to match your local paths
file_zlist = "so2_zlist"
file_ylist = "so2_ylist"
file_xlist = "so2_xlist"

file_vp = "so2.mesh_vp"
file_vs = "so2.mesh_vs"
file_rho = "so2.mesh_rho"

output_nc = "model_CVMSI_taper_so2.nc"

## Read ASCII coordinate files
zlist = np.loadtxt(file_zlist, dtype=np.float32)
ylist = np.loadtxt(file_ylist, dtype=np.float32)
xlist = np.loadtxt(file_xlist, dtype=np.float32)

# Extract dimension sizes dynamically from coordinate lists
n_depth = len(zlist)      # Expecting 210
n_lat = len(ylist)        # Expecting 1251
n_lon = len(xlist)        # Expecting 1301

print(f"n_depth is '{n_depth}'")
print(f"n_lat is '{n_lat}'")
print(f"n_lon is '{n_lon}'")

## Read raw binary 3D data files and reshape to (depth, latitude, longitude)
# Note: Adjust byteorder if needed (e.g., dtype='>f4' for big-endian)
vp_data = np.fromfile(file_vp, dtype=np.float32).reshape((n_depth, n_lat, n_lon))
vs_data = np.fromfile(file_vs, dtype=np.float32).reshape((n_depth, n_lat, n_lon))
rho_data = np.fromfile(file_rho, dtype=np.float32).reshape((n_depth, n_lat, n_lon))
vp_cnt = vp_data.size
vs_cnt = vs_data.size
rho_cnt = rho_data.size

print(f"vp_cnt is '{vp_cnt}'")
print(f"vs_cnt is '{vs_cnt}'")
print(f"rho_cnt is '{rho_cnt}'")

#print(f"vs..'{vs_data[100][100][100]}'")

## Create and write NetCDF4 dataset
with Dataset(output_nc, mode="w", format="NETCDF4") as ncfile:

    # Define dimensions
    ncfile.createDimension("depth", n_depth)
    ncfile.createDimension("latitude", n_lat)
    ncfile.createDimension("longitude", n_lon)

    # Define coordinate variables
    var_depth = ncfile.createVariable("depth", "f4", ("depth",), fill_value=np.nan)
    var_depth.units = "meters"

    var_lat = ncfile.createVariable("latitude", "f4", ("latitude",), fill_value=np.nan)
    var_lat.units = "degrees_north"

    var_lon = ncfile.createVariable("longitude", "f4", ("longitude",), fill_value=np.nan)
    var_lon.units = "degrees_east"

    # Define 3D data variables
    var_vp = ncfile.createVariable("vp", "f4", ("depth", "latitude", "longitude"), fill_value=np.nan)
    var_vp.units = "m/s"

    var_vs = ncfile.createVariable("vs", "f4", ("depth", "latitude", "longitude"), fill_value=np.nan)
    var_vs.units = "m/s"

    var_rho = ncfile.createVariable("rho", "f4", ("depth", "latitude", "longitude"), fill_value=np.nan)
    var_rho.units = "kg/m^3"

    # Assign values
    var_depth[:] = zlist
    var_lat[:] = ylist
    var_lon[:] = xlist

    var_vp[:] = vp_data
    var_vs[:] = vs_data
    var_rho[:] = rho_data

## All done
print(f"NetCDF file '{output_nc}' successfully written.")