Adjust nz calculation to match EMOD3D changes in generate_3d_model.py

[andy22b]

Hey team!

Just tried using this to get VS30 for some sites. I was getting an issue with the depth check in threshold.py, I used Claude to find it and have pasted the issue below. Just using the previous release (v2026.4.1) for now. I also notices there's been a change in syntax since one of the README examples was written.

Thanks heaps!

Nz offbyone issue · MD

generate_1d_profiles and generate_thresholds fail for every station: nz off-by-one vs. gen_full_model_grid_great_circle

Summary

Since e0dbb96 (merged in #59, velocity_model_node_correction), gen_full_model_grid_great_circle expects nz to include the EMOD3D padding row, but generate_1d_profiles.py and threshold.py still compute nz without it. Every station therefore fails validation, independent of coordinates, extent or spacing.

Released version 2026.7.1 on PyPI contains the commit, so a fresh pip install velocity-modelling reproduces this on the first station.

Reproduction

echo "id,lon,lat,zmin,zmax,spacing" > sites.csv
echo "STATION_A,172.724,-43.515,0,2,0.01" >> sites.csv
 
generate_1d_profiles sites.csv --model-version 2.03 --out-dir /tmp/profiles

Actual

RuntimeError: Failed to process station STATION_A at (-43.515, 172.724):  nz: 200 !=  (zmax-zmin: 2.0-0.0) / h_depth: 0.01)(=201)

Expected

A profile is generated.

Cause

velocity_modelling/geometry.py:1095 validates against:

nz_expected = int((zmax - zmin) / h_depth + 0.5) + 1

The + 1 was added in e0dbb96 to match generate_3d_model.py (~L197-203), which adds a padding row because EMOD3D does not read the last layer of the velocity model. That commit touched geometry.py only.

The other two callers that build a mesh through the same function still use the un-padded form:

• velocity_modelling/scripts/generate_1d_profiles.py ~L586:

  model_extent["nz"] = int(
      (model_extent["extent_zmax"] - model_extent["extent_zmin"])
      / model_extent["h_depth"]
      + 0.5
  )

• velocity_modelling/threshold.py:602:

  nz = int((zmax - zmin) / h_depth + 0.5)

So int(2.0 / 0.01 + 0.5) = 200 on the caller side against 201 in the validator. This is not a floating-point edge case — 2.0 / 0.01 is exactly 200.0.

The depth_values branch in generate_1d_profiles.py is affected too: nz is computed the same way with the h_depth = 1.0 placeholder before gen_full_model_grid_great_circle is called, and gm.nz is only overwritten afterwards.

Affected versions

Tag	Contains e0dbb96
v2026.4.1 and earlier	no
v2026.07.1 (current PyPI 2026.7.1)	yes

Suggested fix

The padding row is an EMOD3D output requirement, not a property of the grid geometry, so hard-coding + 1 in the shared validator arguably pushes an output-format concern into geometry.py. Two options:

1. Move the padding back to the caller: revert nz_expected to int((zmax - zmin) / h_depth + 0.5) and keep the + 1 in generate_3d_model.py only, where the EMOD3D rationale already lives.
2. Keep the validator as is and add + 1 in generate_1d_profiles.py and threshold.py.
   Option 2 changes 1D profile output: the profile becomes endpoint-inclusive (0.00-2.00 km rather than 0.00-1.99 km for the example above), which may be desirable but is a behavioural change worth calling out in the changelog. Note also that geometry.py:1130 shifts the shallowest node down by h_depth / 4, so the first sample is at 2.5 m rather than 0 m.

A regression test covering generate_1d_profiles and generate_thresholds end-to-end for a single station would catch this class of drift between the callers and the validator.

[gemini-code-assist]

Caution

The consumer version of Gemini Code Assist on GitHub has been sunset. All code review activity has officially ceased.