diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
new file mode 100644
index 0000000..64937b1
--- /dev/null
+++ b/.github/workflows/pytest.yml
@@ -0,0 +1,24 @@
+name: Pytest Check
+on: [pull_request]
+
+jobs:
+ test:
+ runs-on: ubuntu-latest
+ steps:
+ - name: Remove junk
+ uses: easimon/maximize-build-space@master
+ with:
+ root-reserve-mb: 512
+ swap-size-mb: 1024
+ remove-dotnet: "true"
+ - name: Checkout code
+ uses: nschloe/action-cached-lfs-checkout@v1
+
+ - name: Setup UV
+ uses: astral-sh/setup-uv@v5
+ with:
+ enable-cache: true
+ cache-dependency-glob: "**/pyproject.toml"
+
+ - name: Run Tests
+ run: uv run --extra test pytest tests
diff --git a/.github/workflows/ruff.yml b/.github/workflows/ruff.yml
new file mode 100644
index 0000000..8e10796
--- /dev/null
+++ b/.github/workflows/ruff.yml
@@ -0,0 +1,8 @@
+name: Ruff
+on: [pull_request]
+jobs:
+ ruff:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: chartboost/ruff-action@v1
diff --git a/.github/workflows/types.yml b/.github/workflows/types.yml
new file mode 100644
index 0000000..26351b2
--- /dev/null
+++ b/.github/workflows/types.yml
@@ -0,0 +1,22 @@
+name: Type Check
+
+on: [pull_request]
+
+jobs:
+ typecheck:
+ runs-on: ubuntu-latest
+
+ steps:
+ - name: Checkout code
+ uses: actions/checkout@v4
+
+ - name: Setup Python
+ uses: actions/setup-python@v5
+
+ - name: Install uv
+ uses: astral-sh/setup-uv@v5
+ with:
+ enable-cache: true
+
+ - name: Run type checking with ty
+ run: uv run --all-extras ty check
diff --git a/.github/workflows/yamllint.yml b/.github/workflows/yamllint.yml
new file mode 100644
index 0000000..c35dfdb
--- /dev/null
+++ b/.github/workflows/yamllint.yml
@@ -0,0 +1,17 @@
+name: "Yamllint GitHub Actions"
+on:
+ - pull_request
+jobs:
+ yamllint:
+ name: "Yamllint"
+ runs-on: ubuntu-latest
+ steps:
+ - name: "Checkout"
+ uses: actions/checkout@master
+ - name: "Yamllint"
+ uses: karancode/yamllint-github-action@master
+ with:
+ yamllint_strict: true
+ yamllint_comment: true
+ env:
+ GITHUB_ACCESS_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/.yamllint b/.yamllint
new file mode 100644
index 0000000..c5b0b01
--- /dev/null
+++ b/.yamllint
@@ -0,0 +1,3 @@
+rules:
+ line-length:
+ max: 200
diff --git a/nzcvm_registry.yaml b/nzcvm_registry.yaml
index dea1d99..a50b83b 100644
--- a/nzcvm_registry.yaml
+++ b/nzcvm_registry.yaml
@@ -1,34 +1,247 @@
+---
tomography:
- name: EP2010
- elev: [ 15, 1, -3, -8, -15, -23, -30, -38, -48, -65, -85, -105, -130, -155, -185, -225, -275, -370, -620, -750 ]
- path: tomography/EP2010/ep2010.h5
+ path: tomography/EP2010/EP2010.h5
+ elev:
+ [
+ 15,
+ 1,
+ -3,
+ -8,
+ -15,
+ -23,
+ -30,
+ -38,
+ -48,
+ -65,
+ -85,
+ -105,
+ -130,
+ -155,
+ -185,
+ -225,
+ -275,
+ -370,
+ -620,
+ -750,
+ ]
author: Eberhart-Phillips et al. (2010)
title: Establishing a Versatile 3-d seismic Velocity Model for New Zealand
- url: https://10.1785/gssrl.81.6.992
+ url: https://www.seismosoc.org/Publications/srl/SRL_81/srl_81-6_eberhart-phillips_et_al-esupp/Table_S1.txt
+
+ - name: EP2017
+ path: tomography/EP2017/EP2017.h5
+ elev:
+ [
+ 15,
+ 1,
+ -1,
+ -3,
+ -5,
+ -8,
+ -15,
+ -23,
+ -30,
+ -34,
+ -38,
+ -42,
+ -48,
+ -55,
+ -65,
+ -85,
+ -105,
+ -130,
+ -155,
+ -185,
+ -225,
+ -275,
+ -370,
+ -620,
+ -750,
+ ]
+ author: Eberhart-Phillips et al. (2017)
+ title: New Zealand Wide model 2.1 seismic velocity and Qs and Qp models for New Zealand
+ url: https://zenodo.org/record/1043558
- name: EP2020
- elev: [ 15, 1, -1, -3, -5, -8, -15, -23, -30, -34, -38, -42, -48, -55, -65, -85, -105, -130, -155, -185, -225, -275, -370, -620, -750 ]
- path: tomography/EP2020/ep2020.h5
+ path: tomography/EP2020/EP2020.h5
+ elev:
+ [
+ 15,
+ 1,
+ -1,
+ -3,
+ -5,
+ -8,
+ -15,
+ -23,
+ -30,
+ -34,
+ -38,
+ -42,
+ -48,
+ -55,
+ -65,
+ -85,
+ -105,
+ -130,
+ -155,
+ -185,
+ -225,
+ -275,
+ -370,
+ -620,
+ -750,
+ ]
author: Eberhart-Phillips et al. (2020)
title: New Zealand Wide model 2.2 seismic velocity and Qs and Qp models for New Zealand
- url: https://10.5281/zenodo.3779523
-
-
- - name: CHOW2020_EP2020_MIX
- elev: [ 15.0, 2.25, 2.0, 1.75, 1.5, 1.25, 1.0, 0.75, 0.5, 0.25, 0.0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -1.75, -2.0, -2.25, -2.5, -2.75, -3.0, -3.25, -3.5, -3.75, -4.0, -4.25, -4.5, -4.75, -5.0, -5.25, -5.5, -5.75, -6.0, -6.25, -6.5, -6.75, -7.0, -7.25, -7.5, -7.75, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0, -17.0, -18.0, -19.0, -20.0, -21.0, -22.0, -23.0, -24.0, -25.0, -26.0, -27.0, -28.0, -29.0, -30.0, -31.0, -32.0, -33.0, -34.0, -35.0, -36.0, -37.0, -38.0, -39.0, -40.0, -41.0, -42.0, -43.0, -44.0, -45.0, -46.0, -47.0, -48.0, -49.0, -50.0, -52.0, -56.0, -60.0, -64.0, -68.0, -72.0, -76.0, -80.0, -84.0, -88.0, -92.0, -96.0, -100.0, -104.0, -108.0, -112.0, -116.0, -120.0, -124.0, -128.0, -132.0, -136.0, -140.0, -144.0, -148.0, -152.0, -156.0, -160.0, -164.0, -168.0, -172.0, -176.0, -180.0, -184.0, -188.0, -192.0, -196.0, -200.0, -204.0, -208.0, -212.0, -216.0, -220.0, -224.0, -228.0, -232.0, -236.0, -240.0, -244.0, -248.0, -252.0, -256.0, -260.0, -264.0, -268.0, -272.0, -276.0, -280.0, -284.0, -288.0, -292.0, -296.0, -300.0, -304.0, -308.0, -312.0, -316.0, -320.0, -324.0, -328.0, -332.0, -336.0, -340.0, -344.0, -348.0, -352.0, -356.0, -360.0, -364.0, -368.0, -372.0, -376.0, -380.0, -384.0, -388.0, -392.0, -396.0, -400.0, -620.0, -750.0, ]
- path: tomography/CHOW2020_EP2020_MIX/chow2020_ep2020_mix.h5
- author: Chow et al. (2020)
- url:
- - https://doi.org/10.1093/gji/ggaa381
- - https://core.geo.vuw.ac.nz/d/feae69f61ea54f81bee1
+ url: https://zenodo.org/records/3779523
+
+ - name: EP2022
+ path: tomography/EP2022/EP2022.h5
+ elev:
+ [
+ 15,
+ 1,
+ -1,
+ -3,
+ -5,
+ -8,
+ -15,
+ -23,
+ -30,
+ -34,
+ -38,
+ -42,
+ -48,
+ -55,
+ -65,
+ -85,
+ -105,
+ -130,
+ -155,
+ -185,
+ -225,
+ -275,
+ -370,
+ -620,
+ -750,
+ ]
+ author: Eberhart-Phillips et al. (2022)
+ title: New Zealand Wide model 2.3 seismic velocity and Qs and Qp models for New Zealand
+ url: https://zenodo.org/records/5098356
- name: EP2025
- elev: [ 1, -3, -8, -15, -23, -30, -38, -48, -65, -85, -105, -130, -155, -185, -225, -275, -370, -620, -750 ]
- path: tomography/EP2025/ep2025.h5
+ path: tomography/EP2025/EP2025.h5
+ elev:
+ [
+ 15,
+ 1,
+ -1,
+ -3,
+ -5,
+ -8,
+ -15,
+ -23,
+ -30,
+ -34,
+ -38,
+ -42,
+ -48,
+ -55,
+ -65,
+ -85,
+ -105,
+ -130,
+ -155,
+ -185,
+ -225,
+ -275,
+ -370,
+ -620,
+ -750,
+ ]
author: Eberhart-Phillips et al. (2025)
- title: New Zealand Wide model 2.2 seismic velocity and Qs and Qp models for New Zealand
- url: n/a
+ title: New Zealand Wide model 3.1 seismic velocity and Qs and Qp models for New Zealand
+ url: Personal communication (pending publication)
+
+ - name: SW2025
+ path: tomography/SW2025/SW2025_EP2022_SmoothConvolve.h5
+ elev:
+ [
+ 15, 5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5, -6, -7, -8,
+ -9, -10, -11, -12, -13, -14, -15, -16, -17, -18, -19,
+ -20, -21, -22, -23, -24, -25, -26, -27, -28, -29,
+ -30, -31, -32, -33, -34, -35, -36, -37, -38, -39,
+ -40, -41, -42, -43, -44, -45, -46, -47, -48, -49,
+ -50, -51, -52, -53, -54, -55, -56, -57, -58, -59,
+ -60, -65, -85, -105, -130, -155, -185, -225, -275,
+ -370, -620, -750
+ ]
+ author: Wu et al. (2025)
+ title: Seismic azimuthal anisotropy of New Zealand revealed by adjoint-state traveltime tomography
+ url: Full Model Data - Personal Communication, Downsampled data published in https://doi.org/10.21979/N9/QKU80S
+
+ - name: DB2025
+ path: tomography/DB2025/DB2025_EP2025_SmoothConvolve.h5
+ elev:
+ [
+ 15.0, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5, -0.0, -0.5, -1.0,
+ -1.5, -2.0, -2.5, -3.0, -3.5, -4.0, -4.5, -5.0, -5.5, -6.0,
+ -6.5, -7.0, -7.5, -8.0, -8.5, -9.0, -9.5, -10.0, -10.5, -11.0,
+ -11.5, -12.0, -12.5, -13.0, -13.5, -14.0, -14.5, -15.0, -15.5, -16.0,
+ -16.5, -17.0, -17.5, -18.0, -18.5, -19.0, -19.5, -20.0, -20.5, -21.0,
+ -21.5, -22.0, -22.5, -23.0, -23.5, -24.0, -24.5, -25.0, -25.5, -26.0,
+ -26.5, -27.0, -27.5, -28.0, -28.5, -29.0, -29.5, -30.0, -30.5, -31.0,
+ -31.5, -32.0, -32.5, -33.0, -33.5, -34.0, -34.5, -35.0, -35.5, -36.0,
+ -36.5, -37.0, -37.5, -38.0, -38.5, -39.0, -39.5, -40.0, -40.5, -41.0,
+ -41.5, -42.0, -42.5, -43.0, -43.5, -44.0, -44.5, -45.0, -45.5, -46.0,
+ -46.5, -47.0, -47.5, -48.0, -48.5, -49.0, -49.5, -50.0, -50.5, -51.0,
+ -51.5, -52.0, -52.5, -53.0, -53.5, -54.0, -54.5, -55.0, -55.5, -56.0,
+ -56.5, -57.0, -57.5, -58.0, -58.5, -59.0, -59.5, -60.0, -60.5, -61.0,
+ -61.5, -62.0, -62.5, -63.0, -63.5, -64.0, -64.5, -65.0, -65.5, -66.0,
+ -66.5, -67.0, -67.5, -68.0, -68.5, -69.0, -69.5, -70.0, -70.5, -71.0,
+ -71.5, -72.0, -72.5, -73.0, -73.5, -74.0, -74.5, -75.0, -75.5, -76.0,
+ -76.5, -77.0, -77.5, -78.0, -78.5, -79.0, -79.5, -80.0, -80.5, -81.0,
+ -81.5, -82.0, -82.5, -83.0, -83.5, -84.0, -84.5, -85.0, -85.5, -86.0,
+ -86.5, -87.0, -87.5, -88.0, -88.5, -89.0, -89.5, -90.0, -90.5, -91.0,
+ -91.5, -92.0, -92.5, -93.0, -93.5, -94.0, -94.5, -95.0, -95.5, -96.0,
+ -96.5, -97.0, -97.5, -98.0, -98.5, -99.0, -99.5, -100.0,
+ -105.0, -130.0, -155.0, -185.0, -225.0, -275.0, -370.0, -620.0, -750.0
+ ]
+ author: Bassett et al. (2025)
+ title: Crustal Structure of the Hikurangi Subduction Zone Revealed by Four Decades of Onshore-Offshore Seismic Data
+ url: Latest model data obtained through personal communication (email). Older data published in https://zenodo.org/records/13381669
+
+ - name: BC2022
+ path: tomography/BC2022/BC2022_EP2020_SmoothConvolve.h5
+ elev:
+ [
+ 15.0, 2.25, 2.0, 1.75, 1.5, 1.25, 1.0, 0.75, 0.5, 0.25,
+ 0.0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -1.75, -2.0, -2.25,
+ -2.5, -2.75, -3.0, -3.25, -3.5, -3.75, -4.0, -4.25, -4.5, -4.75,
+ -5.0, -5.25, -5.5, -5.75, -6.0, -6.25, -6.5, -6.75, -7.0, -7.25,
+ -7.5, -7.75, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0,
+ -16.0, -17.0, -18.0, -19.0, -20.0, -21.0, -22.0, -23.0, -24.0, -25.0,
+ -26.0, -27.0, -28.0, -29.0, -30.0, -31.0, -32.0, -33.0, -34.0, -35.0,
+ -36.0, -37.0, -38.0, -39.0, -40.0, -41.0, -42.0, -43.0, -44.0, -45.0,
+ -46.0, -47.0, -48.0, -49.0, -50.0, -52.0, -55.0, -56.0, -60.0, -64.0,
+ -65.0, -68.0, -72.0, -76.0, -80.0, -84.0, -85.0, -88.0, -92.0, -96.0,
+ -100.0, -104.0, -105.0, -108.0, -112.0, -116.0, -120.0, -124.0, -128.0, -130.0,
+ -132.0, -136.0, -140.0, -144.0, -148.0, -152.0, -155.0, -156.0, -160.0, -164.0,
+ -168.0, -172.0, -176.0, -180.0, -184.0, -185.0, -188.0, -192.0, -196.0, -200.0,
+ -204.0, -208.0, -212.0, -216.0, -220.0, -224.0, -225.0, -228.0, -232.0, -236.0,
+ -240.0, -244.0, -248.0, -252.0, -256.0, -260.0, -264.0, -268.0, -272.0, -275.0,
+ -276.0, -280.0, -284.0, -288.0, -292.0, -296.0, -300.0, -304.0, -308.0, -312.0,
+ -316.0, -320.0, -324.0, -328.0, -332.0, -336.0, -340.0, -344.0, -348.0, -352.0,
+ -356.0, -360.0, -364.0, -368.0, -370.0, -372.0, -376.0, -380.0, -384.0, -388.0,
+ -392.0, -396.0, -400.0, -620.0, -750.0
+ ]
+ author: Chow et al. (2022)
+ title: Strong Upper-Plate Heterogeneity at the Hikurangi Subduction Margin (North Island, New Zealand) Imaged by Adjoint Tomography
+ url: Latest model data obtained through https://doi.org/10.17611/dp/emc.2021.nzatomnnorthvpvs.1.
+
basin:
- name: Canterbury_v18p1
boundaries:
@@ -42,8 +255,7 @@ basin:
submodel: miocene_submod_v1
- path: regional/Canterbury/Canterbury_Paleogene_WGS84.h5
submodel: paleogene_submod_v1
- - path: regional/Canterbury/Canterbury_Basement_WGS84.h5
-
+ - path: regional/Canterbury/Canterbury_basement_WGS84.h5
- name: Canterbury_v18p2
boundaries:
@@ -57,8 +269,7 @@ basin:
submodel: miocene_submod_v1
- path: regional/Canterbury/Canterbury_Paleogene_WGS84.h5
submodel: paleogene_submod_v1
- - path: regional/Canterbury/Canterbury_Basement_WGS84.h5
-
+ - path: regional/Canterbury/Canterbury_basement_WGS84.h5
- name: Canterbury_v18p3
boundaries:
@@ -72,8 +283,7 @@ basin:
submodel: miocene_submod_v1
- path: regional/Canterbury/Canterbury_Paleogene_WGS84.h5
submodel: paleogene_submod_v1
- - path: regional/Canterbury/Canterbury_Basement_WGS84.h5
-
+ - path: regional/Canterbury/Canterbury_basement_WGS84.h5
- name: Canterbury_v19p1
type: 4
@@ -121,7 +331,6 @@ basin:
- path: regional/Canterbury/Canterbury_basement_WGS84.h5
smoothing: regional/Canterbury/Canterbury_smoothing.txt
-
- name: NorthCanterbury_v19p1
type: 1
author: Robin Lee
@@ -135,7 +344,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/NorthCanterbury/NorthCanterbury_basement_WGS84_v19p1.h5
-
- name: NorthCanterbury_v25p8
type: 1
author: Robin Lee / Ayushi Tiwari
@@ -164,7 +372,6 @@ basin:
submodel: bpv_submod_v4
- path: regional/BanksPeninsulaVolcanics/BanksPeninsulaVolcanics_Miocene_WGS84.h5
-
- name: Kaikoura_v19p1
type: 2
author: Robin Lee
@@ -180,7 +387,6 @@ basin:
- path: regional/Kaikoura/Kaikoura_basement_WGS84.h5
smoothing: regional/Kaikoura/Kaikoura_smoothing.txt
-
- name: Kaikoura_v25p5
type: 2
author: Robin Lee
@@ -195,7 +401,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Kaikoura/Kaikoura_basement_WGS84_v25p5.h5
smoothing: regional/Kaikoura/Kaikoura_smoothing_v25p5.txt
-
- name: Cheviot_v19p1
type: 1
@@ -211,7 +416,6 @@ basin:
- path: regional/Cheviot/Cheviot_basement_WGS84.h5
smoothing: regional/Cheviot/Cheviot_smoothing.txt
-
- name: Hanmer_v19p1
type: 1
author: Robin Lee
@@ -226,7 +430,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Hanmer/Hanmer_basement_WGS84_v19p1.h5
-
- name: Hanmer_v25p3
type: 1
author: Ayushi Tiwari
@@ -254,7 +457,6 @@ basin:
- path: regional/Marlborough/Marlborough_basement_WGS84.h5
smoothing: regional/Marlborough/Marlborough_smoothing.txt
-
- name: Nelson_v19p1
type: 2
author: Robin Lee
@@ -271,7 +473,6 @@ basin:
- path: regional/Nelson/Nelson_basement_WGS84_v19p1.h5
smoothing: regional/Nelson/Nelson_smoothing.txt
-
- name: Nelson_v25p5
type: 3
author: Robin Lee
@@ -306,7 +507,6 @@ basin:
- path: regional/Nelson/Nelson_basement_WGS84_v25p5.h5
smoothing: regional/Nelson/Nelson_smoothing.txt
-
- name: Wellington_v19p1
boundaries:
- regional/Wellington/Wellington_outline_WGS84.geojson
@@ -316,17 +516,15 @@ basin:
- path: regional/Wellington/Wellington_basement_WGS84_v19p1.h5
smoothing: regional/Wellington/Wellington_smoothing.txt
-
- name: Wellington_v19p6
boundaries:
- regional/Wellington/Wellington_outline_WGS84.geojson
surfaces:
- path: surface/NZ_DEM_HD.h5
submodel: canterbury1d_v2
- - path : regional/Wellington/Wellington_basement_WGS84_v19p6.h5
+ - path: regional/Wellington/Wellington_basement_WGS84_v19p6.h5
smoothing: regional/Wellington/Wellington_smoothing.txt
-
- name: Wellington_v21p8
type: 3
author: Robin Lee / Matt Hill
@@ -379,7 +577,6 @@ basin:
- path: regional/WaikatoHauraki/WaikatoHauraki_basement_WGS84.h5
smoothing: regional/WaikatoHauraki/WaikatoHauraki_smoothing.txt
-
- name: Wanaka_v20p6
type: 1
author: Cameron Douglas (USER2020)
@@ -394,8 +591,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Wanaka/Wanaka_basement_WGS84.h5
-
-
- name: Mackenzie_v20p6
type: 1
author: Cameron Douglas (USER2020)
@@ -410,7 +605,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Mackenzie/Mackenzie_basement_WGS84.h5
-
- name: Wakatipu_v20p7
type: 1
author: Tim Tuckey (USER2020)
@@ -424,7 +618,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Wakatipu/Wakatipu_basement_WGS84.h5
-
- name: Alexandra_v20p7
type: 1
author: Cameron Douglas (USER2020)
@@ -441,7 +634,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Alexandra/Alexandra_basement_WGS84.h5
-
- name: Ranfurly_v20p7
type: 1
author: Cameron Douglas (USER2020)
@@ -458,7 +650,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Ranfurly/Ranfurly_basement_WGS84.h5
-
- name: NE_Otago_v20p7
type: 1
author: Cameron Douglas (USER2020)
@@ -479,7 +670,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/NE_Otago/NE_Otago_basement_WGS84.h5
-
- name: Mosgiel_v20p7
type: 1
author: Cameron Douglas (USER2020)
@@ -500,8 +690,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Mosgiel/Mosgiel_basement_WGS84.h5
-
-
- name: Balclutha_v20p7
type: 1
author: Cameron Douglas (USER2020)
@@ -540,8 +728,6 @@ basin:
- path: regional/Dunedin/Dunedin_basement_WGS84.h5
smoothing: regional/Dunedin/Dunedin_smoothing.txt
-
-
- name: Murchison_v20p7
type: 1
author: Tim Tuckey (USER2020)
@@ -556,8 +742,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Murchison/Murchison_basement_WGS84.h5
-
-
- name: Waitaki_v20p8
type: 1
author: Cameron Douglas (USER2020)
@@ -577,7 +761,6 @@ basin:
- path: regional/Waitaki/Waitaki_basement_WGS84.h5
smoothing: regional/Waitaki/Waitaki_smoothing.txt
-
- name: Hakataramea_v20p8
type: 1
author: Cameron Douglas (USER2020)
@@ -594,7 +777,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/Hakataramea/Hakataramea_basement_WGS84.h5
-
- name: Karamea_v20p11
type: 1
author: Tim Tuckey (USER2020)
@@ -609,7 +791,6 @@ basin:
- path: regional/Karamea/Karamea_basement_WGS84.h5
smoothing: regional/Karamea/Karamea_smoothing.txt
-
- name: Collingwood_v20p11
type: 1
author: Tim Tuckey (USER2020)
@@ -626,7 +807,6 @@ basin:
- path: regional/Collingwood/Collingwood_basement_WGS84.h5
smoothing: regional/Collingwood/Collingwood_smoothing.txt
-
- name: SpringsJunction_v20p11
type: 1
author: Tim Tuckey (USER2020)
@@ -640,7 +820,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/SpringsJunction/SpringsJunction_basement_WGS84.h5
-
- name: HawkesBay_v21p7
type: 1
author: William Lee (USER2021)
@@ -724,7 +903,6 @@ basin:
wiki_images:
- images/NI_mideast.png
- images/Gisborne_basin_map.png
- - images/grisborne_basement.png
notes:
- (Comment from the author)
- Stop at Moto River, working around the East Cape from Gisborne. Haurere Point should be included in Whakatane
@@ -736,7 +914,6 @@ basin:
- path: regional/Gisborne/Gisborne_basement_WGS84.h5
smoothing: regional/Gisborne/Gisborne_smoothing.txt
-
- name: SouthernHawkesBay_v21p12
type: 1
author: William Lee (USER2021)
@@ -757,7 +934,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/SouthernHawkesBay/SouthernHawkesBay_basement_WGS84.h5
-
- name: Wairarapa_v21p12
type: 1
author: William Lee (USER2021)
@@ -765,7 +941,7 @@ basin:
- images/NI_south.png
- images/Wairarapa_basin_map.png
notes:
- - (Comment from the author) "Consider adding east coastal basins (e.g. Uruti Point)"
+ - (Comment from the author) "Consider adding east coastal basins (e.g. Uruti Point)"
boundaries:
- regional/Wairarapa/Wairarapa_outline_WGS84.geojson
surfaces:
@@ -774,7 +950,6 @@ basin:
- path: regional/Wairarapa/Wairarapa_basement_WGS84.h5
smoothing: regional/Wairarapa/Wairarapa_smoothing.txt
-
- name: OmaioBay_v22p3
type: 1
author: Cameron Davis / Emma Coumbe (USER2022)
@@ -798,7 +973,6 @@ basin:
- path: regional/OmaioBay/OmaioBay_basement_WGS84_v22p3.h5
smoothing: regional/OmaioBay/OmaioBay_smoothing_v22p3.txt
-
- name: OmaioBay_v25p5
type: 1
author: Cameron Davis / Emma Coumbe (USER2022)
@@ -821,8 +995,6 @@ basin:
- path: regional/OmaioBay/OmaioBay_basement_WGS84.h5
smoothing: regional/OmaioBay/OmaioBay_smoothing.txt
-
-
- name: Whangaparoa_v23p4
type: 1
author: Cameron Davis / Emma Coumbe (USER2022)
@@ -838,8 +1010,7 @@ basin:
- path: surface/NZ_DEM_HD.h5
submodel: canterbury1d_v2
- path: regional/Whangaparoa/Whangaparoa_basement_WGS84.h5
- smoothing: regional/Whangaparoa/Whangaparoa_smoothing.txt
-
+ smoothing: regional/Whangaparoa/Whangaparoa_smoothing.txt
- name: PalmerstonNorth_v25p5
type: 1
@@ -847,11 +1018,11 @@ basin:
wiki_images:
- images/PalmerstonNorth_basin_map_v25p5.png
boundaries:
- - regional/PalmerstonNorth/PalmerstonNorth_outline_WGS84.geojson
+ - regional/PalmerstonNorth/PalmerstonNorth_outline_WGS84.geojson
surfaces:
- - path: surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path: regional/PalmerstonNorth/PalmerstonNorth_basement_WGS84_v25p5.h5
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/PalmerstonNorth/PalmerstonNorth_basement_WGS84_v25p5.h5
smoothing: regional/PalmerstonNorth/PalmerstonNorth_smoothing.txt
- name: PalmerstonNorth_v25p8
@@ -865,11 +1036,10 @@ basin:
- path: surface/NZ_DEM_HD.h5
submodel: palmerstonnorth_v1
- path: regional/PalmerstonNorth/PalmerstonNorth_pliocenetop_WGS84_v25p8.h5
- submodel: pn_pliocene_submod_v1
+ submodel: pn_pliocene_submod_v1
- path: regional/PalmerstonNorth/PalmerstonNorth_basement_WGS84_v25p8.h5
smoothing: regional/PalmerstonNorth/PalmerstonNorth_smoothing.txt
-
- name: Southland_v25p5
type: 2
author: Archie Goodrick
@@ -898,7 +1068,6 @@ basin:
submodel: canterbury1d_v2
- path: regional/TeAnau/TeAnau_basement_WGS84.h5
-
- name: TolagaBay_v25p5
type: 1
author: Cameron Davis / Emma Coumbe (USER2022)
@@ -908,7 +1077,7 @@ basin:
boundaries:
- regional/TolagaBay/TolagaBay_outline_WGS84.geojson
-
+
surfaces:
- path: surface/NZ_DEM_HD.h5
submodel: canterbury1d_v2
@@ -929,7 +1098,6 @@ basin:
- path: regional/Waiapu/Waiapu_basement_WGS84.h5
smoothing: regional/Waiapu/Waiapu_smoothing.txt
-
- name: WestCoast_v25p5
type: 1
author: Ayushi Tiwari / Hunter Brotherston
@@ -938,12 +1106,11 @@ basin:
boundaries:
- regional/WestCoast/WestCoast_outline_WGS84.geojson
surfaces:
- - path: surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path: regional/WestCoast/WestCoast_basement_WGS84.h5
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/WestCoast/WestCoast_basement_WGS84.h5
smoothing: regional/WestCoast/WestCoast_smoothing.txt
-
- name: Westport_v25p5
type: 1
author: Ayushi Tiwari / Kaleb Finn (ENCN493)
@@ -970,155 +1137,150 @@ basin:
- path: regional/Westport/Westport_basement_WGS84.h5
smoothing: regional/Westport/Westport_smoothing.txt
-
- name: QueenCharlotte_v25p8
type: 1
author: Ayushi Tiwari
wiki_images:
- - images/QueenCharlotte_basin_map.png
+ - images/QueenCharlotte_basin_map.png
boundaries:
- - regional/QueenCharlotte/QueenCharlotte_outline_WGS84.geojson
+ - regional/QueenCharlotte/QueenCharlotte_outline_WGS84.geojson
surfaces:
- - path: surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path: regional/QueenCharlotte/QueenCharlotte_basement_WGS84_v25p8.h5
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/QueenCharlotte/QueenCharlotte_basement_WGS84_v25p8.h5
- name: CastleHill_v25p8
type: 1
author: Ayushi Tiwari
wiki_images:
- - images/CastleHill_basin_map.png
+ - images/CastleHill_basin_map.png
boundaries:
- - regional/CastleHill/CastleHill_outline_WGS84.geojson
+ - regional/CastleHill/CastleHill_outline_WGS84.geojson
surfaces:
- - path: surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path: regional/CastleHill/CastleHill_basement_WGS84.h5
-
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/CastleHill/CastleHill_basement_WGS84.h5
- name: Whakatane_v25p8
type: 1
author: Ayushi Tiwari
wiki_images:
- - images/Whakatane_basin_map.png
+ - images/Whakatane_basin_map.png
boundaries:
- - regional/Whakatane/Whakatane_outline_WGS84.geojson
+ - regional/Whakatane/Whakatane_outline_WGS84.geojson
surfaces:
- - path: surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path: regional/Whakatane/Whakatane_basement_WGS84_v25p8.h5
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/Whakatane/Whakatane_basement_WGS84_v25p8.h5
smoothing: regional/Whakatane/Whakatane_smoothing.txt
- name: TeAraroa_v25p9
type: 1
author: James Hoskin
wiki_images:
- - images/TeAraroa_basin_map.png
+ - images/TeAraroa_basin_map.png
boundaries:
- - regional/TeAraroa/TeAraroa_outline_WGS84.geojson
+ - regional/TeAraroa/TeAraroa_outline_WGS84.geojson
surfaces:
- - path : surface/NZ_DEM_HD.h5
+ - path: surface/NZ_DEM_HD.h5
submodel: canterbury1d_v2
- - path : regional/TeAraroa/TeAraroa_basement_WGS84.h5
+ - path: regional/TeAraroa/TeAraroa_basement_WGS84.h5
smoothing: regional/TeAraroa/TeAraroa_smoothing.txt
- name: Rarakau_v25p9
type: 1
author: Toby Bates
wiki_images:
- - images/Rarakau_basin_map.png
+ - images/Rarakau_basin_map.png
boundaries:
- - regional/Rarakau/Rarakau_outline_WGS84.geojson
+ - regional/Rarakau/Rarakau_outline_WGS84.geojson
surfaces:
- - path : surface/NZ_DEM_HD.h5
- submodel: canterbury1d_v2
- - path : regional/Rarakau/Rarakau_basement_WGS84.h5
+ - path: surface/NZ_DEM_HD.h5
+ submodel: canterbury1d_v2
+ - path: regional/Rarakau/Rarakau_basement_WGS84.h5
smoothing: regional/Rarakau/Rarakau_smoothing.txt
-
vs30:
- name: nz_with_offshore
path: vs30/NZ_Vs30_HD_With_Offshore.h5
- name: nz_without_offshore
path: vs30/NZ_Vs30.h5
-
submodel:
- - name: nan_submod
- type: null
-
- - name: canterbury1d_v1
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/Cant1D_v1.fd_modfile
-
- - name: canterbury1d_v2
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/Cant1D_v2.fd_modfile
-
- - name: canterbury1d_v2_pliocene_enforced
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/Cant1D_v2_Pliocene_Enforced.fd_modfile
-
- - name: canterbury1d_v3_pliocene_enforced
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/Cant1D_v3_Pliocene_Enforced.fd_modfile
-
- - name: nelson_v1
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/Nelson_v1.fd_modfile
-
- - name: palmerstonnorth_v1
- type: vm1d
- module: canterbury1d_submod
- data: vm1d/PalmerstonNorth_v1.fd_modfile
-
- - name: ep_tomography_submod_v2010
- type: tomography
- module: ep_tomography_submod_v2010
-
- - name: bpv_submod_v4
- type: relation
- module: bpv_submod_v4
-
- - name: pliocene_submod_v1
- type: relation
- module: pliocene_submod_v1
-
- - name: pn_pliocene_submod_v1
- type: relation
- module: pn_pliocene_submod_v1
-
-
- - name: miocene_submod_v1
- type: relation
- module: miocene_submod_v1
-
- - name: paleogene_submod_v1
- type: relation
- module: paleogene_submod_v1
-
- - name: pliocene_submod_v2
- type: relation
- module: pliocene_submod_v2
-
- - name: miocene_submod_v2
- type: relation
- module: miocene_submod_v2
-
- - name: paleogene_submod_v2
- type: relation
- module: paleogene_submod_v2
-
- - name: perturbation_v20p6
- type: perturbation
-
- - name: perturbation_v20p10
- type: perturbation
-
- - name: perturbation_v20p11
- type: perturbation
+ - name: nan_submod
+ type: null
+
+ - name: canterbury1d_v1
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/Cant1D_v1.fd_modfile
+
+ - name: canterbury1d_v2
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/Cant1D_v2.fd_modfile
+
+ - name: canterbury1d_v2_pliocene_enforced
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/Cant1D_v2_Pliocene_Enforced.fd_modfile
+
+ - name: canterbury1d_v3_pliocene_enforced
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/Cant1D_v3_Pliocene_Enforced.fd_modfile
+
+ - name: nelson_v1
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/Nelson_v1.fd_modfile
+
+ - name: palmerstonnorth_v1
+ type: vm1d
+ module: canterbury1d_submod
+ data: vm1d/PalmerstonNorth_v1.fd_modfile
+
+ - name: ep_tomography_submod_v2010
+ type: tomography
+ module: ep_tomography_submod_v2010
+
+ - name: bpv_submod_v4
+ type: relation
+ module: bpv_submod_v4
+
+ - name: pliocene_submod_v1
+ type: relation
+ module: pliocene_submod_v1
+
+ - name: pn_pliocene_submod_v1
+ type: relation
+ module: pn_pliocene_submod_v1
+
+ - name: miocene_submod_v1
+ type: relation
+ module: miocene_submod_v1
+
+ - name: paleogene_submod_v1
+ type: relation
+ module: paleogene_submod_v1
+
+ - name: pliocene_submod_v2
+ type: relation
+ module: pliocene_submod_v2
+
+ - name: miocene_submod_v2
+ type: relation
+ module: miocene_submod_v2
+
+ - name: paleogene_submod_v2
+ type: relation
+ module: paleogene_submod_v2
+
+ - name: perturbation_v20p6
+ type: perturbation
+
+ - name: perturbation_v20p10
+ type: perturbation
+
+ - name: perturbation_v20p11
+ type: perturbation
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..91ed9f8
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,21 @@
+
+[project]
+name = "nzcvm_data"
+requires-python = ">=3.11"
+version = "0.1.0"
+description = "Integrity checks for NZCVM data registry"
+dependencies = ["numpy>=2.3.5", "pytest-subtests>=0.15.0", "schema>=0.7.8"]
+
+
+[project.optional-dependencies]
+test = [
+ "h5py>=3.15.1",
+ "pytest>=9.0.2",
+ "pyyaml>=6.0.3",
+ "geojson>=3.2.0",
+ "shapely>=2.1.2",
+]
+# Even though this isn't a python package, these dev dependencies can
+# still prove useful for writing tests.
+dev = ["ruff", "numpydoc", "ty"]
+types = ["types-shapely>=2.1.0.20250917"]
diff --git a/regional/Dunedin/Dunedin_smoothing.txt b/regional/Dunedin/Dunedin_smoothing.txt
index 9499c95..9faae13 100755
--- a/regional/Dunedin/Dunedin_smoothing.txt
+++ b/regional/Dunedin/Dunedin_smoothing.txt
@@ -131,21 +131,3 @@
170.45759988 -45.92441961
170.45719142 -45.92356637
170.45689042 -45.92269490
-170.56320143 -45.86537325
-170.56356933 -45.86624843
-170.56393725 -45.86712360
-170.56430518 -45.86799877
-170.56467313 -45.86887395
-170.56504108 -45.86974911
-170.56540905 -45.87062428
-170.56577703 -45.87149945
-170.56614502 -45.87237461
-170.56651302 -45.87324977
-170.56688103 -45.87412493
-170.56724906 -45.87500009
-170.56761709 -45.87587525
-170.56798514 -45.87675041
-170.56835320 -45.87762556
-170.56872127 -45.87850071
-170.56908935 -45.87937587
-170.56945745 -45.88025101
diff --git a/tests/test_registry.py b/tests/test_registry.py
new file mode 100644
index 0000000..a5fe59a
--- /dev/null
+++ b/tests/test_registry.py
@@ -0,0 +1,584 @@
+import itertools
+import re
+from pathlib import Path
+from typing import NotRequired, TypedDict, no_type_check
+
+import h5py
+import numpy as np
+import pytest
+import shapely
+import yaml
+from attr import dataclass
+from pytest import Metafunc
+from pytest_subtests import SubTests
+from schema import Optional, Or, Regex, Schema, SchemaError
+
+
+@pytest.fixture(scope="session")
+def nzcvm_registry_path() -> Path:
+ return Path(__file__).parent.parent / "nzcvm_registry.yaml"
+
+
+@pytest.fixture(scope="session")
+def nzcvm_root() -> Path:
+ return Path(__file__).parent.parent
+
+
+@no_type_check
+def test_nzcvm_registry_schema(nzcvm_registry_path: Path) -> None:
+ path = Regex(
+ # See https://stackoverflow.com/a/537876
+ r"[^\0]+",
+ error="Must be valid unix path.",
+ )
+ ident = Regex(r"^[a-zA-Z_][a-zA-Z0-9_]*$", error="Must be valid python identifier.")
+ # Source - https://stackoverflow.com/a
+ # Posted by Daveo, modified by community. See post 'Timeline' for change history
+ # Retrieved 2025-12-17, License - CC BY-SA 4.0
+ url = Or(
+ Regex(
+ r"https?://(www\.)?[-a-zA-Z0-9@:%._\+~#=]{1,256}\.[a-zA-Z0-9()]{1,6}\b([-a-zA-Z0-9()@:%_\+.~#?&//=]*)"
+ ),
+ "Personal communication (pending publication)",
+ error='Must be a valid URL, or "Personal communication (pending publication)"',
+ )
+ surface_schema = Schema({"path": path, Optional("submodel"): ident})
+
+ tomography_entry = Schema(
+ {
+ "name": ident,
+ "elev": [Or(int, float)],
+ "path": path,
+ "author": str,
+ Optional("title"): str,
+ Optional("url"): Or(
+ url,
+ [url],
+ ), # Handles single URL, list of URLs, or empty
+ }
+ )
+
+ basin_entry = Schema(
+ {
+ "name": ident,
+ Optional("type"): Or(1, 2, 3, 4),
+ Optional("author"): str,
+ Optional("notes"): [str],
+ Optional("wiki_images"): [path],
+ "boundaries": [path],
+ "surfaces": [surface_schema],
+ Optional("smoothing"): str,
+ }
+ )
+
+ submodel_schema = Schema(
+ {
+ "name": ident,
+ "type": Or("vm1d", "tomography", "relation", "perturbation", None),
+ Optional("module"): ident,
+ Optional("data"): path,
+ }
+ )
+ vs30_schema = Schema({"name": str, "path": path})
+ registry_schema = Schema(
+ {
+ "tomography": [tomography_entry],
+ "basin": [basin_entry],
+ "submodel": [submodel_schema],
+ "vs30": [vs30_schema],
+ }
+ )
+
+ with open(nzcvm_registry_path, "r") as f:
+ registry = yaml.safe_load(f)
+
+ try:
+ registry_schema.validate(registry)
+ except SchemaError as e:
+ pytest.fail(f"NZCVM Registry is invalid\n{e.code}")
+
+
+class Tomography(TypedDict):
+ name: str
+ elev: list[float]
+ path: str
+ author: str
+ title: str
+ url: str
+
+
+class Surface(TypedDict):
+ path: str
+ submodel: NotRequired[str]
+
+
+class Basin(TypedDict):
+ name: str
+ author: str
+ notes: list[str]
+ wiki_images: list[str]
+ boundaries: list[str]
+ surfaces: list[Surface]
+ smoothing: NotRequired[str]
+
+
+class Vs30(TypedDict):
+ name: str
+ path: str
+
+
+class Submodel(TypedDict):
+ name: str
+ type: str
+ module: str
+ data: NotRequired[str]
+
+
+def pytest_generate_tests(metafunc: Metafunc) -> None:
+ registry_path = Path(__file__).parent.parent / "nzcvm_registry.yaml"
+ with open(registry_path) as f:
+ registry = yaml.safe_load(f)
+ if "model" in metafunc.fixturenames:
+ # This assumes you have access to the registry here
+ # It creates a separate test case for every model entry
+ metafunc.parametrize("model", registry["tomography"], ids=lambda m: m["name"])
+ elif "basin" in metafunc.fixturenames:
+ metafunc.parametrize("basin", registry["basin"], ids=lambda m: m["name"])
+ elif "vs30" in metafunc.fixturenames:
+ metafunc.parametrize("vs30", registry["vs30"], ids=lambda m: m["name"])
+ elif "submodel" in metafunc.fixturenames:
+ metafunc.parametrize("submodel", registry["submodel"], ids=lambda m: m["name"])
+
+
+def test_tomography_paths_exist(nzcvm_root: Path, model: Tomography) -> None:
+ relative_model_path = Path(model["path"])
+ model_path = nzcvm_root / relative_model_path
+ assert model_path.exists()
+
+
+def test_tomography_models_are_hdf5(nzcvm_root: Path, model: Tomography) -> None:
+ relative_model_path = Path(model["path"])
+ model_path = nzcvm_root / relative_model_path
+ try:
+ with h5py.File(model_path, "r") as f:
+ assert len(f.keys()) > 0
+ except Exception as e:
+ pytest.fail(f"{model['name']} is not a valid hdf5 file: {e}")
+
+
+def test_tomography_elevations_match(nzcvm_root: Path, model: Tomography) -> None:
+ relative_model_path = Path(model["path"])
+ model_path = nzcvm_root / relative_model_path
+ registry_elevations = sorted(model["elev"])
+ with h5py.File(model_path, "r") as f:
+ model_elevations = sorted([float(elev) for elev in f.keys()])
+
+ assert registry_elevations == model_elevations, "Model elevations do not match."
+
+
+def test_tomography_compatible_shapes(
+ subtests: SubTests, nzcvm_root: Path, model: Tomography
+) -> None:
+ model_path = nzcvm_root / model["path"]
+
+ with h5py.File(model_path, "r") as f:
+ for elev, group in f.items():
+ name = model["name"]
+ with subtests.test(msg=f"Checking elevation {elev}", model=name, elev=elev):
+ lat_shape = group["latitudes"].shape
+ lon_shape = group["longitudes"].shape
+
+ # Check data arrays against (lat, lon)
+ expected_data_shape = (lat_shape[0], lon_shape[0])
+
+ for field in ["vp", "vs", "rho"]:
+ actual_shape = group[field].shape
+ assert actual_shape == expected_data_shape, (
+ f"Shape mismatch in {model_path.name} at {elev}m: "
+ f"{field} is {actual_shape}, expected {expected_data_shape}"
+ )
+
+
+R_EARTH = 6378.139 # km, from: https://github.com/ucgmsim/velocity_modelling/blob/27c7e6e64d7ce1a9e543d58a6e584d498358431c/velocity_modelling/constants.py#L15
+LONGITUDE_TOLERANCE = 0.01 # km
+LATITUDE_TOLERANCE = 0.01 # km
+LAT_DEGREES_PER_KM = np.pi / 180 * R_EARTH
+
+
+def test_tomography_geo_gridpoints(
+ subtests: SubTests, nzcvm_root: Path, model: Tomography
+) -> None:
+ relative_model_path = Path(model["path"])
+ model_path = nzcvm_root / relative_model_path
+
+ with h5py.File(model_path, "r") as f:
+ for elev, group in f.items():
+ name = model["name"]
+ with subtests.test(msg=f"Checking elevation {elev}", model=name, elev=elev):
+ latitude = np.array(group["latitudes"])
+ longitude = np.array(group["longitudes"])
+
+ lat_diffs_km = np.diff(latitude) * LAT_DEGREES_PER_KM
+ assert np.all(lat_diffs_km > 0), "Latitudes not strictly ascending"
+ assert latitude[0] >= -90 and latitude[-1] <= 90, (
+ "Latitudes must be between -90 and 90."
+ )
+ assert lat_diffs_km == pytest.approx(
+ np.full(len(lat_diffs_km), lat_diffs_km[0]), abs=LATITUDE_TOLERANCE
+ )
+
+ lon_diffs_deg = np.diff(longitude) # Shape (N-1,)
+ assert np.all(lon_diffs_deg > 0), "Longitudes not strictly ascending"
+ assert longitude[0] >= 0 and longitude[-1] <= 185, (
+ "Longitudes must be between 0 and 185."
+ )
+
+
+# NOTE: Values for QUALITY_BOUNDS are not physically derived
+QUALITY_BOUNDS = {"vp": (0, 11.0), "vs": (0, 7.0), "rho": (0, 5.0)}
+
+
+@pytest.mark.parametrize("quality", ["vp", "rho", "vs"])
+def test_tomography_quality(
+ subtests: SubTests, nzcvm_root: Path, model: Tomography, quality: str
+) -> None:
+ relative_model_path = Path(model["path"])
+ model_path = nzcvm_root / relative_model_path
+
+ with h5py.File(model_path, "r") as f:
+ for elev, group in f.items():
+ name = model["name"]
+ with subtests.test(msg=f"Checking elevation {elev}", model=name, elev=elev):
+ quality_values = np.array(group[quality])
+ assert not np.isnan(quality_values).any(), (
+ f"Quality {quality} contains NaN values."
+ )
+ bounds = QUALITY_BOUNDS.get(quality)
+ assert bounds
+ min, max = bounds
+ quality_min = quality_values.min()
+ quality_max = quality_values.max()
+ assert quality_min >= min, (
+ f"Quality {quality} minimum value ({quality_min=}) is less than {min}."
+ )
+ assert quality_max <= max, (
+ f"Quality {quality} maximum value ({quality_max=}) is less than {max}."
+ )
+
+
+def test_basin_boundaries_exist(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for boundary in basin["boundaries"]:
+ boundary_relative_path = Path(boundary)
+ name = basin["name"]
+ with subtests.test(
+ msg=f"Checking boundary {boundary_relative_path.stem}",
+ basin=name,
+ boundary=boundary_relative_path.stem,
+ ):
+ boundary_path = nzcvm_root / boundary_relative_path
+ assert boundary_path.exists()
+
+
+def test_basin_boundaries_are_valid_geojson(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for boundary in basin["boundaries"]:
+ boundary_relative_path = Path(boundary)
+ name = basin["name"]
+ boundary_name = boundary_relative_path.stem
+
+ with subtests.test(
+ msg=f"Checking boundary {boundary_name}", basin=name, boundary=boundary_name
+ ):
+ boundary_path = nzcvm_root / boundary_relative_path
+
+ try:
+ geojson_str = boundary_path.read_text()
+ geom_collection = shapely.from_geojson(geojson_str)
+ except (OSError, ValueError, Exception) as e:
+ pytest.fail(f"Model {name} has invalid boundary {boundary_name}: {e}")
+
+ assert isinstance(geom_collection, shapely.GeometryCollection)
+ assert shapely.is_valid(geom_collection), (
+ "Geometry is topologically invalid"
+ )
+
+ assert not shapely.is_empty(geom_collection), "Geometry is empty"
+
+ assert all(
+ isinstance(geom, shapely.Polygon) for geom in geom_collection.geoms
+ )
+
+
+def test_basin_paths_exist(nzcvm_root: Path, basin: Basin) -> None:
+ if match := re.match(r"^(\w+)_v\d+p\d+$", basin["name"]):
+ basin_canonical_name = match.group(1)
+ assert basin_canonical_name
+ basin_path = nzcvm_root / "regional" / basin_canonical_name
+ assert basin_path.exists()
+ if "wiki_images" in basin:
+ assert all(
+ (basin_path / Path(image)).exists() for image in basin["wiki_images"]
+ )
+ else:
+ pytest.fail("basin name does not follow structure name_vxxpxx")
+
+
+def test_basin_surfaces_exist(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for surface in basin["surfaces"]:
+ relative_surface_path = Path(surface["path"])
+ name = basin["name"]
+ surface_name = relative_surface_path.stem
+ with subtests.test(
+ msg=f"Checking surface {surface_name}", basin=name, surface=surface_name
+ ):
+ surface_path = nzcvm_root / relative_surface_path
+ assert surface_path.exists()
+
+
+def test_basin_surfaces_are_valid_hdf5(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for surface in basin["surfaces"]:
+ relative_surface_path = Path(surface["path"])
+ name = basin["name"]
+ surface_name = relative_surface_path.stem
+ with subtests.test(
+ msg=f"Checking surface {surface_name}", basin=name, surface=surface_name
+ ):
+ surface_path = nzcvm_root / relative_surface_path
+ try:
+ with h5py.File(surface_path, "r") as f:
+ assert "elevation" in f.keys()
+ assert "latitude" in f.keys()
+ assert "longitude" in f.keys()
+ except Exception as e:
+ pytest.fail(f"{surface_name} is not a valid hdf5 file: {e}")
+
+
+def test_surface_geo_gridpoints(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for surface in basin["surfaces"]:
+ relative_surface_path = Path(surface["path"])
+ name = basin["name"]
+ surface_name = relative_surface_path.stem
+ with subtests.test(
+ msg=f"Checking surface {surface_name}", basin=name, surface=surface_name
+ ):
+ surface_path = nzcvm_root / relative_surface_path
+ with h5py.File(surface_path, "r") as f:
+ latitude = np.array(f["latitude"])
+ longitude = np.array(f["longitude"])
+
+ lat_diffs_km = np.diff(latitude) * LAT_DEGREES_PER_KM
+ assert np.all(lat_diffs_km > 0), "Latitudes not strictly ascending"
+ assert latitude[0] >= -90 and latitude[-1] <= 90, (
+ "Latitudes must be between -90 and 90."
+ )
+
+ lon_diffs_deg = np.diff(longitude) # Shape (N-1,)
+ assert np.all(lon_diffs_deg > 0), "Longitudes not strictly ascending"
+ assert longitude[0] >= 0 and longitude[-1] <= 185, (
+ "Longitudes must be between 0 and 185."
+ )
+
+ elevation = np.array(f["elevation"])
+ assert elevation.shape == (len(latitude), len(longitude)), (
+ "Elevation shape must be match latitude and longitude"
+ )
+ assert not np.isnan(elevation).any(), "Elevations cannot be NaN"
+ assert elevation.min() >= -10000, "Elevations cannot be below -10000m"
+ assert elevation.max() <= 10000, "Elevations cannot be above 10000m"
+
+
+def read_smoothing_boundary(smoothing_path: Path) -> shapely.LineString:
+ coords: list[tuple[float, float]] = []
+ with open(smoothing_path, "r") as f:
+ for line in f:
+ line = line.strip()
+ if not line:
+ continue
+ line_coords = re.split(r"\s+", line)
+ assert len(line_coords) == 2
+ lon_str, lat_str = line_coords
+ assert lon_str and lat_str
+ coords.append((float(lon_str), float(lat_str)))
+
+ return shapely.LineString(coords)
+
+
+def test_basin_smoothing_contained_in_boundaries(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ if "smoothing" not in basin:
+ pytest.skip("basin has no smoothing boundary")
+
+ boundaries = []
+ for boundary in basin["boundaries"]:
+ boundary_relative_path = Path(boundary)
+ boundary_path = nzcvm_root / boundary_relative_path
+ geojson_str = boundary_path.read_text()
+ geom_collection = shapely.from_geojson(geojson_str)
+ assert isinstance(geom_collection, shapely.GeometryCollection)
+ for geom in geom_collection.geoms:
+ assert isinstance(geom, shapely.Polygon)
+ boundaries.append(geom.exterior)
+
+ # Add a small smoothing boundary buffer to account for the fact
+ # that smoothing boundary is not perfectly contained in the basin
+ # boundary.
+ boundary_geometry = shapely.buffer(shapely.union_all(boundaries), 0.001)
+ smoothing_surface_relative_path = basin["smoothing"]
+ smoothing_surface_path = nzcvm_root / Path(smoothing_surface_relative_path)
+ smoothing_boundary = read_smoothing_boundary(smoothing_surface_path)
+
+ assert boundary_geometry.contains(smoothing_boundary)
+
+
+def test_basin_surfaces_contain_boundaries(
+ subtests: SubTests, nzcvm_root: Path, basin: Basin
+) -> None:
+ for boundary, surface in itertools.product(basin["boundaries"], basin["surfaces"]):
+ boundary_relative_path = Path(boundary)
+ surface_relative_path = Path(surface["path"])
+ basin_name = basin["name"]
+ boundary_name = boundary_relative_path.stem
+ surface_name = surface_relative_path.stem
+
+ with subtests.test(
+ msg=f"Checking boundary {boundary_name}",
+ basin=basin_name,
+ boundary=boundary_name,
+ surface=surface_name,
+ ):
+ boundary_path = nzcvm_root / boundary_relative_path
+
+ geojson_str = boundary_path.read_text()
+ geom_collection = shapely.from_geojson(geojson_str)
+ surface_path = nzcvm_root / surface_relative_path
+ with h5py.File(surface_path, "r") as f:
+ latitudes = np.array(f["latitude"])
+ longitudes = np.array(f["longitude"])
+ elevation_boundary = shapely.box(
+ xmin=longitudes.min(),
+ xmax=longitudes.max(),
+ ymin=latitudes.min(),
+ ymax=latitudes.max(),
+ )
+ assert shapely.contains(elevation_boundary, geom_collection)
+
+
+def test_vs30_file_exists(nzcvm_root: Path, vs30: Vs30) -> None:
+ relative_path = Path(vs30["path"])
+
+ vs30_path = nzcvm_root / relative_path
+ assert vs30_path.exists(), f"Vs30 file missing at {vs30_path}"
+
+
+def test_vs30_is_valid_hdf5(nzcvm_root: Path, vs30: Vs30) -> None:
+ relative_path = Path(vs30["path"])
+
+ vs30_path = nzcvm_root / relative_path
+
+ try:
+ with h5py.File(vs30_path, "r") as f:
+ assert "elevation" in f.keys(), (
+ "Dataset 'elevation' (vs30) missing from HDF5"
+ )
+ assert "latitude" in f.keys()
+ assert "longitude" in f.keys()
+ except Exception as e:
+ name = vs30["name"]
+ pytest.fail(f"Vs30 file {name} is not a valid hdf5 file: {e}")
+
+
+def test_vs30_geo_gridpoints(nzcvm_root: Path, vs30: Vs30) -> None:
+ relative_path = Path(vs30["path"])
+ vs30_path = nzcvm_root / relative_path
+
+ with h5py.File(vs30_path, "r") as f:
+ latitude = np.array(f["latitude"])
+ longitude = np.array(f["longitude"])
+ vs30_values = np.array(f["elevation"])
+
+ lat_diffs = np.diff(latitude)
+ assert np.all(lat_diffs > 0), "Latitudes not strictly ascending"
+ assert latitude[0] >= -90 and latitude[-1] <= 90, (
+ "Latitudes out of world bounds"
+ )
+
+ lon_diffs = np.diff(longitude)
+ assert np.all(lon_diffs > 0), "Longitudes not strictly ascending"
+ assert longitude[0] >= 0 and longitude[-1] <= 185, (
+ "Longitudes out of NZCVM bounds"
+ )
+
+ assert vs30_values.shape == (len(latitude), len(longitude)), (
+ f"Vs30 shape {vs30_values.shape} does not match lat/lon dimensions"
+ )
+
+ assert not np.isnan(vs30_values).any(), "Vs30 data contains NaNs"
+ assert vs30_values.min() >= 0, (
+ f"Vs30 values below 0 detected: {vs30_values.min()}"
+ )
+ assert vs30_values.max() <= 2000, (
+ f"Vs30 values above 2000 detected: {vs30_values.max()}"
+ )
+
+
+def test_submodel_data_exists_where_relevant(
+ nzcvm_root: Path, submodel: Submodel
+) -> None:
+ assert ("data" in submodel) == (submodel["type"] == "vm1d"), (
+ "Submodel has data iff submodel is a vm1d"
+ )
+ if "data" in submodel:
+ data_relative_path = Path(submodel["data"])
+ data_path = nzcvm_root / data_relative_path
+ assert data_path.exists()
+
+
+@dataclass
+class SubmodelData:
+ vp: float
+ vs: float
+ rho: float
+ qp: float
+ qs: float
+ thickness: float
+
+
+def parse_submodel_data(submodel_path: Path) -> list[SubmodelData]:
+ rows = []
+ with open(submodel_path, "r") as f:
+ header = next(f).strip()
+ assert header == "DEF HST"
+ for line in f:
+ row = re.split(r"\s+", line.strip())
+ floats = [float(x) for x in row]
+ assert len(floats) == 6
+ rows.append(SubmodelData(*floats))
+ return rows
+
+
+def test_submodel_data_is_valid(
+ subtests: SubTests, nzcvm_root: Path, submodel: Submodel
+) -> None:
+ if "data" in submodel:
+ data_relative_path = Path(submodel["data"])
+ data_path = nzcvm_root / data_relative_path
+ data = parse_submodel_data(data_path)
+ vp_min, vp_max = QUALITY_BOUNDS["vp"]
+ assert all(vp_min <= row.vp <= vp_max for row in data)
+ vs_min, vs_max = QUALITY_BOUNDS["vs"]
+ assert all(vs_min <= row.vs <= vs_max for row in data)
+ assert all(row.thickness >= 0 for row in data)
+ assert all(row.qp > 0 for row in data)
+ assert all(row.qs > 0 for row in data)
+ else:
+ pytest.skip(f"Submodel {submodel['name']} has no data")
diff --git a/tomography/BASSETT2025/README.md b/tomography/BASSETT2025/README.md
deleted file mode 100644
index 117a71e..0000000
--- a/tomography/BASSETT2025/README.md
+++ /dev/null
@@ -1,5 +0,0 @@
-This is currently in progress.
-
-----
-Bassett, D., Henrys, S., Tozer, B., van Avendonk, H., Gase, A., Bangs, N., et al. (2025). Crustal structure of the Hikurangi subduction zone revealed by four decades of Onshore-Offshore seismic data: Implications for the dimensions and slip behavior of the seismogenic zone. Journal of Geophysical Research: Solid Earth, 130, e2024JB030268. https://doi.org/10.1029/2024JB030268
-
diff --git a/tomography/BC2022/BC2022_EP2020_SmoothConvolve.h5 b/tomography/BC2022/BC2022_EP2020_SmoothConvolve.h5
new file mode 100644
index 0000000..8814b08
--- /dev/null
+++ b/tomography/BC2022/BC2022_EP2020_SmoothConvolve.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:135af5a0b9486417032eeab170592956e35877a96876e1eaec4603e3813dba47
+size 3004682172
diff --git a/tomography/CHOW2020_EP2020_MIX/chow2020_ep2020_mix.h5 b/tomography/CHOW2020_EP2020_MIX/chow2020_ep2020_mix.h5
deleted file mode 100644
index a82f872..0000000
--- a/tomography/CHOW2020_EP2020_MIX/chow2020_ep2020_mix.h5
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:bc5600a90df9806662a3afa2f6e109b93de0ac76df5dc08512bc20618a0eb51b
-size 734537257
diff --git a/tomography/DB2025/DB2025_EP2025_SmoothConvolve.h5 b/tomography/DB2025/DB2025_EP2025_SmoothConvolve.h5
new file mode 100644
index 0000000..457f29b
--- /dev/null
+++ b/tomography/DB2025/DB2025_EP2025_SmoothConvolve.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ee162d6b075d9bdceaa334d52ff89e7bfeacc1aef3d788860aaf485c2233ee4d
+size 3488267598
diff --git a/tomography/EP2010/EP2010.h5 b/tomography/EP2010/EP2010.h5
new file mode 100644
index 0000000..b5835ad
--- /dev/null
+++ b/tomography/EP2010/EP2010.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f16409bb32d5ddbfbaa1e67b7718468acaa589d75a95830585f5143578e2ac1d
+size 498347472
diff --git a/tomography/EP2010/ep2010.h5 b/tomography/EP2010/ep2010.h5
deleted file mode 100644
index 4a2b450..0000000
--- a/tomography/EP2010/ep2010.h5
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:f83d94050e192d7516a9be4cb177303053748f6fe0195147001c128712dc440b
-size 69006576
diff --git a/tomography/EP2017/EP2017.h5 b/tomography/EP2017/EP2017.h5
new file mode 100644
index 0000000..780dd4e
--- /dev/null
+++ b/tomography/EP2017/EP2017.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:60a9d8558b8e59383229bc87ffabc7ce89563dc5f157cc791b99551fefd2fb81
+size 686263471
diff --git a/tomography/EP2020/EP2020.h5 b/tomography/EP2020/EP2020.h5
new file mode 100644
index 0000000..97abb1b
--- /dev/null
+++ b/tomography/EP2020/EP2020.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d22aa31d0b1100c821262e1565d989ea885f4829d03bb440de28893e2d0408d9
+size 691979068
diff --git a/tomography/EP2020/README.md b/tomography/EP2020/README.md
deleted file mode 100644
index 02ccaaa..0000000
--- a/tomography/EP2020/README.md
+++ /dev/null
@@ -1,29 +0,0 @@
-# EP2020 Tomography Model
-
-The EP2020 tomography model, as defined in the NZCVM registry, provides a foundational background velocity structure for New Zealand. It is an essential component for integrating more detailed regional and basin models.
-
-## Model Details
-
-This model is explicitly defined in the `nzcvm_registry.yaml` file, which acts as the central source of truth for its properties:
-
-* **Name:** `EP2020`
-* **Author:** Eberhart-Phillips et al. (2020)
-* **Title:** New Zealand Wide model 2.2 seismic velocity and Q structure
-* **Data Path:** `global/tomography/ep2020.h5`
-* **Elevation Layers (`elev`):** The model is defined on 25 discrete elevation layers, specified in kilometers. These layers range from 15 km above sea level down to 750 km below, providing a deep velocity profile.
-
-## Data Integration and Visualization
-
-The raw data for the EP2020 model consists of scattered points of seismic velocity measurements. For use within the 3D velocity model, this data is processed by interpolating these points onto a regular, uniform grid. The interpolated data forms a series of smooth velocity planes at different elevations.
-
-| Original Data (TXT) | Interpolated Data (HDF5) |
-|---------------------|--------------------------|
-| 
| 
|
-
-The left panel shows the original tomography dataset in EP2020 ASCII format. Grid points follow the model's rotated coordinate system, producing an irregular pattern of longitudes that cross the dateline.
-
-The right panel shows the same model after interpolation onto a uniform rectilinear latitude–longitude grid. The interpolated dataset contains 1,120,000 points arranged as 1,400 unique latitudes × 800 unique longitudes, spanning 48°S to 33°S in latitude and 165°E to 180°E in longitude. This regularized grid provides a contiguous New Zealand domain that is directly compatible with 3D visualization and numerical simulations.
-
-## References
-
-* Eberhart-Phillips, D., & Reyners, M. et al. (2020). *New Zealand Wide model 2.2 seismic velocity and Q structure*. New Zealand Journal of Geology and Geophysics.
diff --git a/tomography/EP2020/ep2020.h5 b/tomography/EP2020/ep2020.h5
deleted file mode 100644
index 901da47..0000000
--- a/tomography/EP2020/ep2020.h5
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:62f6fe774a467d5e38bb09d88ed8e42afd86070f32d98d748d523bf660705c7a
-size 110617994
diff --git a/tomography/EP2020/images/ep2020_interpolated_spatial_distribution.png b/tomography/EP2020/images/ep2020_interpolated_spatial_distribution.png
deleted file mode 100644
index 8812927..0000000
Binary files a/tomography/EP2020/images/ep2020_interpolated_spatial_distribution.png and /dev/null differ
diff --git a/tomography/EP2020/images/ep2020_original_spatial_distribution.png b/tomography/EP2020/images/ep2020_original_spatial_distribution.png
deleted file mode 100644
index 148f232..0000000
Binary files a/tomography/EP2020/images/ep2020_original_spatial_distribution.png and /dev/null differ
diff --git a/tomography/EP2020/tools/interpolate_to_uniform_nzcvm_grid.py b/tomography/EP2020/tools/interpolate_to_uniform_nzcvm_grid.py
deleted file mode 100644
index 6f28ee2..0000000
--- a/tomography/EP2020/tools/interpolate_to_uniform_nzcvm_grid.py
+++ /dev/null
@@ -1,200 +0,0 @@
-"""
-Interpolate EP2020 tomography data onto a uniform NZCVM grid.
-
-This script interpolates scattered EP2020 tomography points onto a uniform
-rectilinear grid suitable for use in the NZCVM framework.
-"""
-
-from pathlib import Path
-
-import h5py
-import numpy as np
-import pandas as pd
-from scipy.interpolate import griddata
-
-
-def read_ep_txt(txt_path: str | Path) -> pd.DataFrame:
- """
- Read EP-style TXT tomography data.
-
- Parameters
- ----------
- txt_path : str or Path
- Path to the EP TXT file.
-
- Returns
- -------
- pd.DataFrame
- DataFrame containing tomography data with depth inverted to negative values.
- """
- col_names = [
- "vp", "vp_o_vs", "vs", "rho", "sf_vp", "sf_vp_o_vs",
- "x", "y", "depth", "lat", "lon"
- ]
- df = pd.read_csv(
- txt_path,
- sep=r"\s+",
- skiprows=2,
- names=col_names,
- engine="python"
- )
- df['depth']=-1*df['depth']
-
- return df
-
-
-def load_nzcvm_grid(h5_path: str | Path) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
- """
- Load NZCVM grid coordinates from reference HDF5 file.
-
- Parameters
- ----------
- h5_path : str or Path
- Path to the reference HDF5 file.
-
- Returns
- -------
- tuple[np.ndarray, np.ndarray, np.ndarray]
- Tuple containing (latitudes, longitudes, depths).
- """
- with h5py.File(h5_path, "r") as f:
- depth_keys = sorted(f.keys(), key=lambda x: float(x))
- depth = np.array([float(k) for k in depth_keys])
- group0 = f[depth_keys[0]]
- lat = group0["latitudes"][:]
- lon = group0["longitudes"][:]
- return lat, lon, depth
-
-
-def interpolate_property_to_grid(
- df: pd.DataFrame,
- lat: np.ndarray,
- lon: np.ndarray,
- depth: np.ndarray,
- field_name: str
-) -> np.ndarray:
- """
- Interpolate a property from scattered points to a regular grid.
-
- Parameters
- ----------
- df : pd.DataFrame
- Input DataFrame containing scattered data points.
- lat : np.ndarray
- Target latitude grid points.
- lon : np.ndarray
- Target longitude grid points.
- depth : np.ndarray
- Target depth levels.
- field_name : str
- Name of the field to interpolate.
-
- Returns
- -------
- np.ndarray
- Interpolated data on regular grid with shape (ndepth, nlat, nlon).
- """
- nlat, nlon, ndepth = len(lat), len(lon), len(depth)
- out = np.full((ndepth, nlat, nlon), np.nan, dtype=np.float32)
-
- for iz, d in enumerate(depth):
- df_d = df[np.isclose(df["depth"], d, atol=1e-3)]
- if df_d.empty:
- print(f"⚠️ No data at depth={d:.2f} km")
- continue
-
- points = df_d[["lon", "lat"]].values
- values = df_d[field_name].values
- lon_grid, lat_grid = np.meshgrid(lon, lat)
-
- interp = griddata(points, values, (lon_grid, lat_grid), method="linear", fill_value=0.0)
- out[iz] = interp
-
- return out
-
-
-def write_epstyle_hdf5(
- output_path: str | Path,
- lats: np.ndarray,
- lons: np.ndarray,
- depth_list: np.ndarray,
- vp_stack: np.ndarray,
- vs_stack: np.ndarray,
- rho_stack: np.ndarray
-) -> None:
- """
- Write EP-style tomography HDF5 file.
-
- Each depth slice becomes a group named by depth in km.
-
- Parameters
- ----------
- output_path : str or Path
- Output HDF5 file path.
- lats : np.ndarray
- Latitude values (nlat,).
- lons : np.ndarray
- Longitude values (nlon,).
- depth_list : np.ndarray
- Depths in km (nz,).
- vp_stack : np.ndarray
- P-wave velocity data (nz, nlat, nlon).
- vs_stack : np.ndarray
- S-wave velocity data (nz, nlat, nlon).
- rho_stack : np.ndarray
- Density data (nz, nlat, nlon).
- """
- lon_grid, lat_grid = np.meshgrid(lons, lats)
- coords = np.stack([lat_grid, lon_grid], axis=-1) # shape: (nlat, nlon, 2)
-
- vp_stack[np.isnan(vp_stack)]=-999.0
- vs_stack[np.isnan(vs_stack)]=-999.0
- rho_stack[np.isnan(rho_stack)]=-999.0
-
- with h5py.File(output_path, "w") as f:
- for iz, depth in enumerate(depth_list):
- grp = f.create_group(f"{depth:.0f}")
- grp.create_dataset("latitudes", data=lats)
- grp.create_dataset("longitudes", data=lons)
- grp.create_dataset("coords", data=coords)
- grp.create_dataset("vp", data=vp_stack[iz])
- grp.create_dataset("vs", data=vs_stack[iz])
- grp.create_dataset("rho", data=rho_stack[iz])
-
- print(f"✅ Saved EP-style HDF5 to {output_path}")
-
-
-def main() -> None:
- """
- Main function to run the interpolation process.
-
- This function orchestrates the entire interpolation workflow:
- 1. Read EP2020 data
- 2. Load NZCVM grid
- 3. Interpolate data to grid
- 4. Save results to HDF5
- """
- input_txt = Path("vlnzw2p2dnxyzltln.tbl.txt")
- ref_h5 = Path("../../EP2010/ep2010.h5")
- out_h5 = Path("../ep2020_uniform.h5")
-
- print("📥 Reading EP-style TXT...")
- df = read_ep_txt(input_txt)
- print(set(df['depth']))
-
- print("📐 Loading NZCVM grid...")
- lat, lon, depth = load_nzcvm_grid(ref_h5)
-
- print("📊 Interpolating vp...")
- vp = interpolate_property_to_grid(df, lat, lon, depth, "vp")
- print("📊 Interpolating vs...")
- vs = interpolate_property_to_grid(df, lat, lon, depth, "vs")
- print("📊 Interpolating rho...")
- rho = interpolate_property_to_grid(df, lat, lon, depth, "rho")
-
- print("💾 Writing to output HDF5...")
- write_epstyle_hdf5(out_h5, lat, lon, depth, vp, vs, rho)
-
-
-if __name__ == "__main__":
- main()
diff --git a/tomography/EP2022/EP2022.h5 b/tomography/EP2022/EP2022.h5
new file mode 100644
index 0000000..379fbeb
--- /dev/null
+++ b/tomography/EP2022/EP2022.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f3fb2d96fe6dfcf889bf9a0c63fa265f66b72a0795652a979d1ae2792b7619b5
+size 692517900
diff --git a/tomography/EP2025/EP2025.h5 b/tomography/EP2025/EP2025.h5
new file mode 100644
index 0000000..f18f410
--- /dev/null
+++ b/tomography/EP2025/EP2025.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:79cacdd2bd1f5cb04a65af4b4f8d5ee35e79b31d46adf0e33bbe993b752d00b3
+size 643563147
diff --git a/tomography/EP2025/ep2025.h5 b/tomography/EP2025/ep2025.h5
deleted file mode 100644
index 656e498..0000000
--- a/tomography/EP2025/ep2025.h5
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:82d566ac728581fdad87590d669898621f61f7ccdc079d0432af7e8e65225246
-size 98114796
diff --git a/tomography/README.md b/tomography/README.md
deleted file mode 100644
index 2804d55..0000000
--- a/tomography/README.md
+++ /dev/null
@@ -1,60 +0,0 @@
-# Tomography Models in NZCVM
-
-The tomography models in the NZCVM provide the background velocity structure for New Zealand. These models are derived from seismic travel-time data and offer a lower-resolution (~10 km) representation of the subsurface velocity structure.
-
-## Available Tomography Models
-
-The NZCVM currently supports the following tomography models:
-
-1. EP2010: Based on New Zealand-wide model 1.0 by Eberhart-Phillips et al. (2010).
-2. EP2020: Based on NZWide 2.2 model by Eberhart-Phillips et al. (2020).
-3. EP2025: Based on NZWide 3.1 model by Eberhart-Phillips et al. (2025).
-4. CHOW2020_EP2020_MIX: Combination of the CHOW2020 model (Chow et al. 2020) in North Island and the EP2020 model for the rest
-
-## Tomography Model Definition
-
-Tomography models are defined in the `nzcvm_registry.yaml` file. Here's an example of a tomography model definition:
-
-```yaml
-tomography:
- - name: EP2020
- elev: [ 15, 1, -3, -8, -15, -23, -30, -38, -48, -65, -85, -105, -130, -155, -185, -225, -275, -370, -620, -750 ]
- path: global/tomography/ep2020.h5
- author: Eberhart-Phillips et al. (2020)
- title: New Zealand Wide model 2.2 seismic velocity and Qs and Qp models for New Zealand
- url: https://10.5281/zenodo.3779523
-```
-
-The key components of a tomography model definition are:
-
-- **name**: Identifier for the tomography model.
-- **elev**: Array of elevation values (in kilometers) for the model.
-- **path**: Path to the tomography model data file (in HDF5 format).
-- **author**: Name of the author(s) or organization that created the model.
-- **title**: Title or description of the model.
-- **url**: URL where the model can be accessed or downloaded.
-
-
-## Tomography Submodels
-
-Tomography submodels are used to compute velocity values at specific locations based on the tomography model data. These submodels are defined in the `nzcvm_registry.yaml` file and are associated with surfaces in the model version configuration.
-
-```yaml
-submodel:
- - name: ep_tomography_submod_v2010
- type: tomography
- module: ep_tomography_submod_v2010
-```
-
-The `module` specifies the name of the accompanying Python code that prescribes how to calculate velocity at locations within the region below the `surface`.
-
-
-## Data Format
-
-Tomography model data is stored in HDF5 format. For details on the structure and contents of these files, see the [Data Formats](DataFormats.md) page.
-
-## References
-
-[//]: # (- Donna Eberhart-Phillips, Martin Reyners, Stephen Bannister, Mark Chadwick, Susan Ellis; Establishing a Versatile 3-D Seismic Velocity Model for New Zealand. *Seismological Research Letters* 2010; 81 (6): 992–1000. doi: [https://doi.org/10.1785/gssrl.81.6.992](https://doi.org/10.1785/gssrl.81.6.992).)
-- Donna Eberhart-Phillips, Stephen Bannister, Martin Reyners, and Stuart Henrys. "New Zealand Wide Model 2.2 Seismic Velocity and Qs and Qp Models for New Zealand". *Zenodo*, May 1, 2020. [https://doi.org/10.5281/zenodo.3779523](https://doi.org/10.5281/zenodo.3779523).
-- Bryant Chow, Yoshihiro Kaneko, Carl Tape, Ryan Modrak, John Townend, An automated workflow for adjoint tomography—waveform misfits and synthetic inversions for the North Island, New Zealand, Geophysical Journal International, Volume 223, Issue 3, December 2020, Pages 1461–1480, https://doi.org/10.1093/gji/ggaa381
diff --git a/tomography/SW2025/SW2025_EP2022_SmoothConvolve.h5 b/tomography/SW2025/SW2025_EP2022_SmoothConvolve.h5
new file mode 100644
index 0000000..8722bd4
--- /dev/null
+++ b/tomography/SW2025/SW2025_EP2022_SmoothConvolve.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:69176adbcaf959a648a97b978d957aa95330e28cc76b9189a19ad5e6ebf5a39c
+size 2628504332
diff --git a/tomography/tools/compare_tomo_h5.py b/tomography/tools/compare_tomo_h5.py
deleted file mode 100644
index 77fec5c..0000000
--- a/tomography/tools/compare_tomo_h5.py
+++ /dev/null
@@ -1,372 +0,0 @@
-#!/usr/bin/env python3
-"""
-Compare two HDF5 tomography files using statistical sampling.
-
-This script performs detailed comparison of HDF5 tomography files by sampling
-points from each depth level and computing various error metrics including
-MAE, relative MAE, RMSE, and maximum absolute difference.
-"""
-
-import sys
-import traceback
-from pathlib import Path
-from typing import Annotated
-
-import h5py
-import numpy as np
-import typer
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-VARS_DEFAULT = ("vp", "vs", "rho")
-
-
-def sort_key_level(s: str):
- """
- Sort key function for depth level strings.
-
- Parameters
- ----------
- s : str
- Level string to convert for sorting.
-
- Returns
- -------
- float or str
- Numeric value if convertible, otherwise original string.
- """
- try:
- return float(s)
- except ValueError:
- return s
-
-
-def list_root_groups(h5: h5py.File) -> list[str]:
- """
- List all root groups in an HDF5 file.
-
- Parameters
- ----------
- h5 : h5py.File
- Open HDF5 file handle.
-
- Returns
- -------
- list[str]
- List of root group names.
- """
- return [k for k in h5.keys() if isinstance(h5[k], h5py.Group)]
-
-
-def summarize_file(path: str) -> list[str]:
- """
- Summarize the structure of an HDF5 file.
-
- Parameters
- ----------
- path : str
- Path to the HDF5 file.
-
- Returns
- -------
- list[str]
- List of root group names in the file.
- """
- with h5py.File(path, "r") as f:
- grps = list_root_groups(f)
- print(f"- {path}: {len(grps)} root groups")
- preview = ", ".join(sorted(grps, key=sort_key_level)[:8])
- if len(grps) > 8:
- preview += ", …"
- print(f" groups: {preview}")
- # peek first group for shapes/dtypes
- if grps:
- gname = sorted(grps, key=sort_key_level)[0]
- g = f[gname]
- for ds in ("vp", "vs", "rho", "latitudes", "longitudes"):
- if ds in g:
- print(f" {gname}/{ds}: shape={g[ds].shape}, dtype={g[ds].dtype}")
- return grps
-
-
-def build_row_plan(
- ny: int, nx: int, n_samples: int, rng: np.random.Generator, replace: bool
-) -> tuple[dict[int, np.ndarray], int]:
- """
- Build a sampling plan for selecting points from a 2D grid.
-
- Parameters
- ----------
- ny : int
- Number of rows in the grid.
- nx : int
- Number of columns in the grid.
- n_samples : int
- Target number of samples to select.
- rng : np.random.Generator
- Random number generator.
- replace : bool
- Whether to sample with replacement.
-
- Returns
- -------
- tuple[dict[int, np.ndarray], int]
- A mapping from row indices to column indices, and total unique points count.
- """
- total = ny * nx
- n = min(n_samples, total)
- lin = rng.choice(total, size=n, replace=replace)
- lin = np.unique(lin) # drop duplicates to keep point selection valid
- jj, ii = np.divmod(lin, nx)
-
- plan: dict[int, list] = {}
- for r, c in zip(jj, ii):
- plan.setdefault(int(r), []).append(int(c))
-
- # sort & unique columns per row
- unique_points = 0
- plan_sorted: dict[int, np.ndarray] = {}
- for r, cols in plan.items():
- arr = np.array(cols, dtype=np.int64)
- arr = np.unique(arr) # increasing order
- plan_sorted[r] = arr
- unique_points += arr.size
- return plan_sorted, unique_points
-
-
-def compare_levels(
- f1: h5py.File,
- f2: h5py.File,
- levels: list[str],
- vars_to_check: list[str],
- n_samples: int,
- rng: np.random.Generator,
- with_replacement: bool,
- eps: float,
-):
- """
- Compare data between two HDF5 files across multiple depth levels.
-
- Parameters
- ----------
- f1 : h5py.File
- First HDF5 file handle.
- f2 : h5py.File
- Second HDF5 file handle.
- levels : list[str]
- List of depth levels to compare.
- vars_to_check : list[str]
- List of variable names to compare.
- n_samples : int
- Number of sample points to use for comparison.
- rng : np.random.Generator
- Random number generator.
- with_replacement : bool
- Whether to sample with replacement.
- eps : float
- Epsilon value for relative error calculations.
- """
- totals_abs = {v: 0.0 for v in vars_to_check}
- totals_rel = {v: 0.0 for v in vars_to_check}
- totals_rmse_sq = {v: 0.0 for v in vars_to_check}
- totals_cnt = {v: 0 for v in vars_to_check}
- totals_max = {v: 0.0 for v in vars_to_check}
-
- for lvl in levels:
- print(f"\nLevel {lvl}:")
- g1, g2 = f1[lvl], f2[lvl]
- ref = next((v for v in vars_to_check if v in g1 and v in g2), None)
- if ref is None:
- print(" [skip] none of requested datasets present in both files here")
- continue
- if g1[ref].shape != g2[ref].shape:
- print(
- f" [skip] shape mismatch for {ref}: {g1[ref].shape} vs {g2[ref].shape}"
- )
- continue
-
- ny, nx = g1[ref].shape
- row_plan, n_eff = build_row_plan(ny, nx, n_samples, rng, with_replacement)
- print(f" sampling {n_eff} unique points across {len(row_plan)} rows")
-
- for name in vars_to_check:
- if name not in g1 or name not in g2:
- print(f" {name}: [missing in one file] skip")
- continue
-
- # stream metrics row-by-row to respect monotonic index rule
- abs_sum = 0.0
- rel_sum = 0.0
- sq_sum = 0.0
- n_count = 0
- max_abs = 0.0
-
- ds1, ds2 = g1[name], g2[name]
- for r, cols in row_plan.items():
- a = ds1[r, cols].astype(np.float64, copy=False)
- b = ds2[r, cols].astype(np.float64, copy=False)
- diff = a - b
- absdiff = np.abs(diff)
- denom = np.maximum(eps, np.abs(a))
-
- abs_sum += float(absdiff.sum())
- rel_sum += float((absdiff / denom).sum())
- sq_sum += float((diff * diff).sum())
- n_count += int(a.size)
- if a.size:
- max_abs = max(max_abs, float(absdiff.max()))
-
- if n_count == 0:
- print(f" {name}: no comparable samples")
- continue
-
- mae = abs_sum / n_count
- relmae = rel_sum / n_count
- rmse = np.sqrt(sq_sum / n_count)
- print(
- f" {name}: MAE={mae:.3e} relMAE={relmae:.3e} RMSE={rmse:.3e} max|Δ|={max_abs:.3e}"
- )
-
- totals_abs[name] += abs_sum
- totals_rel[name] += rel_sum
- totals_rmse_sq[name] += sq_sum
- totals_cnt[name] += n_count
- totals_max[name] = max(totals_max[name], max_abs)
-
- print("\n=== Overall (sample-weighted across processed levels) ===")
- for name in vars_to_check:
- if totals_cnt[name] == 0:
- print(f" {name}: no comparable samples")
- continue
- overall_mae = totals_abs[name] / totals_cnt[name]
- overall_rel = totals_rel[name] / totals_cnt[name]
- overall_rmse = np.sqrt(totals_rmse_sq[name] / totals_cnt[name])
- print(
- f" {name}: MAE={overall_mae:.3e} relMAE={overall_rel:.3e} RMSE={overall_rmse:.3e} max|Δ|={totals_max[name]:.3e}"
- )
-
-
-@app.command()
-def main(
- file_a: Annotated[
- Path,
- typer.Argument(exists=True, dir_okay=False, help="First HDF5 file to compare"),
- ],
- file_b: Annotated[
- Path,
- typer.Argument(exists=True, dir_okay=False, help="Second HDF5 file to compare"),
- ],
- vars: Annotated[
- list[str],
- typer.Option(
- "-v",
- "--vars",
- help=f"Datasets to compare (default: {', '.join(VARS_DEFAULT)})",
- ),
- ] = list(VARS_DEFAULT),
- samples: Annotated[
- int,
- typer.Option(
- "-n",
- "--samples",
- min=1,
- help="Target random points per level (unique after de-dup)",
- ),
- ] = 200_000,
- seed: Annotated[int, typer.Option(help="RNG seed")] = 0,
- first_level_only: Annotated[
- bool,
- typer.Option(
- "--first-level-only", help="Compare only the shallowest common level"
- ),
- ] = False,
- with_replacement: Annotated[
- bool,
- typer.Option(
- "--with-replacement", help="Sample with replacement before unique()"
- ),
- ] = False,
- eps: Annotated[
- float,
- typer.Option(min=0.0, help="Epsilon for relative error denominator"),
- ] = 1e-9,
- debug: Annotated[
- bool,
- typer.Option("--debug", help="Show tracebacks on errors"),
- ] = False,
-) -> None:
- """
- Compare two HDF5 tomography files (row-grouped random sampling per level).
-
- Parameters
- ----------
- file_a : Path
- First HDF5 file to compare.
- file_b : Path
- Second HDF5 file to compare.
- vars : list[str], optional
- Datasets to compare.
- samples : int, optional
- Target random points per level.
- seed : int, optional
- RNG seed.
- first_level_only : bool, optional
- Compare only the shallowest common level.
- with_replacement : bool, optional
- Sample with replacement before unique().
- eps : float, optional
- Epsilon for relative error denominator.
- debug : bool, optional
- Show tracebacks on errors.
- """
- rng = np.random.default_rng(seed)
-
- print("=== Structure check ===")
- try:
- grps_a = summarize_file(str(file_a))
- grps_b = summarize_file(str(file_b))
- except (OSError, IOError, KeyError) as e:
- if debug:
- traceback.print_exc()
- else:
- print(f"Error reading HDF5 files: {e}", file=sys.stderr)
- raise typer.Exit(2)
-
- common = sorted(set(grps_a) & set(grps_b), key=sort_key_level)
- only_a = sorted(set(grps_a) - set(grps_b), key=sort_key_level)
- only_b = sorted(set(grps_b) - set(grps_a), key=sort_key_level)
-
- print("\n=== Level set comparison ===")
- print(f"Common levels: {len(common)}")
- if only_a:
- print(
- f"Only in A ({len(only_a)}): {only_a[:8]}{' …' if len(only_a) > 8 else ''}"
- )
- if only_b:
- print(
- f"Only in B ({len(only_b)}): {only_b[:8]}{' …' if len(only_b) > 8 else ''}"
- )
- if not common:
- print("No common levels — nothing to compare.")
- raise typer.Exit(1)
-
- levels = common[:1] if first_level_only else common
- print(
- f"\nWill compare {len(levels)} level(s): {levels[:8]}{' …' if len(levels) > 8 else ''}"
- )
-
- with h5py.File(str(file_a), "r") as f1, h5py.File(str(file_b), "r") as f2:
- compare_levels(
- f1,
- f2,
- levels,
- vars,
- samples,
- rng,
- with_replacement,
- eps,
- )
-
-
-if __name__ == "__main__":
- app()
diff --git a/tomography/tools/tomo_analysis.py b/tomography/tools/tomo_analysis.py
deleted file mode 100644
index 4ff8abf..0000000
--- a/tomography/tools/tomo_analysis.py
+++ /dev/null
@@ -1,422 +0,0 @@
-"""
-Analyze and plot tomography data from TXT or HDF5 files.
-
-This script provides functionality to analyze tomography data, generate spatial
-distribution plots, check grid consistency, and create various visualizations.
-"""
-
-from pathlib import Path
-from typing import Annotated
-
-import cartopy.crs as ccrs
-import cartopy.feature as cfeature
-import h5py
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-import seaborn as sns
-import typer
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-
-def load_txt_ep_format(txt_file: str) -> pd.DataFrame:
- """
- Load EP-style TXT tomography grid. Keeps longitudes as provided
- (may include values >180 or <0).
-
- Parameters
- ----------
- txt_file : str
- Path to the EP-style TXT tomography file.
-
- Returns
- -------
- pd.DataFrame
- DataFrame containing tomography data with columns for coordinates and properties.
- """
- col_names = [
- "vp", "vp_o_vs", "vs", "rho", "sf_vp", "sf_vp_o_vs",
- "x", "y", "depth", "lat", "lon"
- ]
- df = pd.read_csv(
- txt_file,
- sep=r"\s+",
- skiprows=2,
- names=col_names,
- engine="python"
- )
- # NOTE: do NOT force-convert lon here; plotting handles 0–360 or -180..180
- return df
-
-
-def load_hdf5_data(h5_file: str) -> pd.DataFrame:
- """
- Load a single depth slice from the HDF5 tomography file into a flat DataFrame.
-
- Parameters
- ----------
- h5_file : str
- Path to the HDF5 tomography file.
-
- Returns
- -------
- pd.DataFrame
- DataFrame containing flattened tomography data from the shallowest depth slice.
- """
- with h5py.File(h5_file, "r") as f:
- depth_keys = sorted(f.keys(), key=lambda x: float(x))
- group = f[depth_keys[0]]
-
- lat = group["latitudes"][:]
- lon = group["longitudes"][:]
- vp = group["vp"][:] # shape (nlat, nlon)
- lon_grid, lat_grid = np.meshgrid(lon, lat)
-
- df = pd.DataFrame({
- "lat": lat_grid.ravel().astype(np.float32),
- "lon": lon_grid.ravel().astype(np.float32),
- "vp": vp.ravel().astype(np.float32),
- "depth": float(depth_keys[0])
- })
- return df
-
-
-def check_latlon_consistency(h5_file: str) -> None:
- """
- Verify that all depth groups share identical latitude/longitude arrays.
-
- Parameters
- ----------
- h5_file : str
- Path to the HDF5 tomography file to check.
- """
- with h5py.File(h5_file, "r") as f:
- first_group = next(iter(f.keys()))
- lat_ref = f[first_group]["latitudes"][:]
- lon_ref = f[first_group]["longitudes"][:]
-
- all_good = True
- for depth in f.keys():
- lat = f[depth]["latitudes"][:]
- lon = f[depth]["longitudes"][:]
-
- if not (np.allclose(lat, lat_ref) and np.allclose(lon, lon_ref)):
- print(f"Mismatch in lat/lon at depth: {depth}")
- all_good = False
-
- if all_good:
- print("✅ All latitudes and longitudes are consistent across depths.")
- else:
- print("❌ Inconsistencies found in lat/lon arrays.")
-
-
-def check_grid_spacing(lat: np.ndarray, lon: np.ndarray, tol: float = 1e-6) -> None:
- """
- Report unique spacings for lat/lon and whether they are regular.
- Prints approximate km for spacings.
-
- Parameters
- ----------
- lat : np.ndarray
- Array of latitude values in degrees.
- lon : np.ndarray
- Array of longitude values in degrees.
- tol : float, optional
- Tolerance for determining regularity. Default is 1e-6.
- """
- lat_spacing = np.diff(lat)
- lon_spacing = np.diff(lon)
-
- unique_lat_spacings = np.unique(np.round(lat_spacing, decimals=5))
- unique_lon_spacings = np.unique(np.round(lon_spacing, decimals=5))
-
- is_lat_regular = len(unique_lat_spacings) == 1
- is_lon_regular = len(unique_lon_spacings) == 1
-
- # Latitude conversion (constant)
- lat_km_spacing = np.round(unique_lat_spacings * 111.32, 3)
-
- # Longitude spacing in km varies with latitude
- latitudes_for_conversion = np.linspace(lat.min(), lat.max(), num=1000)
- lon_km_spacings = unique_lon_spacings[:, np.newaxis] * 111.32 * np.cos(
- np.radians(latitudes_for_conversion)
- )
- lon_km_min = np.min(lon_km_spacings, axis=1)
- lon_km_max = np.max(lon_km_spacings, axis=1)
-
- print(
- f"Latitude spacings (unique): {unique_lat_spacings} deg ≈ {lat_km_spacing} km | "
- f"Regular: {is_lat_regular}"
- )
- print(
- f"Longitude spacings (unique): {unique_lon_spacings} deg ≈ "
- f"[{lon_km_min.min():.3f} – {lon_km_max.max():.3f}] km | Regular: {is_lon_regular}"
- )
-
-
-def lons_centered_for_display(lon: np.ndarray, center: float = 180.0) -> np.ndarray:
- """
- Shift longitudes into a continuous range [center-180, center+180)
- so there's a single seam at 'center' (default 180°).
-
- Parameters
- ----------
- lon : np.ndarray
- Array of longitude values.
- center : float, optional
- Center longitude for the display range. Default is 180.0.
-
- Returns
- -------
- np.ndarray
- Shifted longitude values.
- """
- lon = np.asarray(lon, dtype=float)
- # wrap to [-180, 180) after subtracting the center, then shift back
- return ((lon - center + 180.0) % 360.0) - 180.0 + center
-
-
-def lons_centered_180(lon: np.ndarray) -> np.ndarray:
- """
- Convert arbitrary longitudes to the native range of
- PlateCarree(central_longitude=180), i.e. [-180, 180).
-
- Parameters
- ----------
- lon : np.ndarray
- Array of longitude values.
-
- Returns
- -------
- np.ndarray
- Converted longitude values in range [-180, 180).
- """
- lon = np.asarray(lon, dtype=float)
- # shift by -180, wrap to [-180,180), done
- return ((lon - 180.0 + 180.0) % 360.0) - 180.0
-
-
-def choose_projection_and_extent(lats: np.ndarray, lons: np.ndarray):
- """
- Decide projection and compute a tight extent.
-
- Parameters
- ----------
- lats : np.ndarray
- Array of latitude values in degrees.
- lons : np.ndarray
- Array of longitude values in degrees.
-
- Returns
- -------
- ax_crs : cartopy.crs.CRS
- Axes projection for plotting.
- data_crs : cartopy.crs.CRS
- CRS of input data (always standard PlateCarree lon/lat).
- extent : tuple
- Extent as (minlon, maxlon, minlat, maxlat) in the CRS given by extent_crs.
- extent_crs : cartopy.crs.CRS
- CRS that 'extent' is expressed in.
- """
- lats = np.asarray(lats, dtype=float)
- lons = np.asarray(lons, dtype=float)
-
- data_crs = ccrs.PlateCarree() # input data CRS
-
- pad_lon = 0.5
- pad_lat = 0.5
- minlat = max(-90.0, float(lats.min()) - pad_lat)
- maxlat = min( 90.0, float(lats.max()) + pad_lat)
-
- has_over_180 = np.nanmax(lons) > 180.0
- has_negative = np.nanmin(lons) < 0.0
- crosses_dl = has_over_180 or has_negative
-
- if crosses_dl:
- # Use 180-centered axes; compute extent in that axes' native range [-180,180)
- ax_crs = ccrs.PlateCarree(central_longitude=180)
- lons_c = lons_centered_180(lons) # <<< key fix (NOT 0..360)
- minlon = max(-180.0, float(lons_c.min()) - pad_lon)
- maxlon = min( 180.0, float(lons_c.max()) + pad_lon)
- extent = (minlon, maxlon, minlat, maxlat)
- extent_crs = ax_crs # extent expressed in axes CRS
- else:
- # Standard case: compute extent in data CRS and pass with data_crs
- ax_crs = ccrs.PlateCarree()
- minlon = max(-180.0, float(lons.min()) - pad_lon)
- maxlon = min( 180.0, float(lons.max()) + pad_lon)
- extent = (minlon, maxlon, minlat, maxlat)
- extent_crs = data_crs
-
- return ax_crs, data_crs, extent, extent_crs
-
-
-def plot_spatial_distribution(df: pd.DataFrame, title_suffix: str = "") -> None:
- """
- Map the grid points without clipping across the dateline.
- Uses a 180°-centered axes when needed and sets extent in the matching CRS.
-
- Parameters
- ----------
- df : pd.DataFrame
- DataFrame containing latitude and longitude columns.
- title_suffix : str, optional
- Suffix to append to the plot title. Default is empty string.
- """
- lats = df["lat"].values
- lons = df["lon"].values
-
- ax_crs, data_crs, extent, extent_crs = choose_projection_and_extent(lats, lons)
-
- plt.figure(figsize=(10, 8))
- ax = plt.axes(projection=ax_crs)
-
- # Set extent in the CRS it was computed in
- ax.set_extent(extent, crs=extent_crs)
-
- ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
- ax.add_feature(cfeature.BORDERS, linestyle=":", linewidth=0.6)
- ax.add_feature(cfeature.LAND, facecolor="lightgray", edgecolor="none")
- ax.add_feature(cfeature.OCEAN, facecolor="lightsteelblue", edgecolor="none")
- ax.add_feature(cfeature.LAKES, alpha=0.5)
- ax.add_feature(cfeature.RIVERS, linewidth=0.5)
-
- ax.scatter(
- df["lon"], df["lat"],
- s=1, c="red", alpha=0.5,
- transform=data_crs, # input data CRS is always standard lon/lat
- label="Model Grid"
- )
- gl = ax.gridlines(draw_labels=True, linestyle="--", linewidth=0.5)
- try:
- gl.right_labels = False
- gl.top_labels = False
- except AttributeError:
- pass
-
- plt.title(f"Spatial Domain of Tomography Model {title_suffix}")
- plt.legend()
- plt.tight_layout()
- plt.show()
-
-
-def plot_spacing_histograms(df: pd.DataFrame) -> None:
- """
- Plot histograms of latitude and longitude spacing and density map.
-
- Parameters
- ----------
- df : pd.DataFrame
- DataFrame containing lat/lon data.
- """
- lat_spacing = np.diff(np.sort(np.unique(df["lat"])))
- lon_spacing = np.diff(np.sort(np.unique(df["lon"])))
-
- # Histogram: latitude spacing
- plt.figure(figsize=(10, 4))
- sns.histplot(lat_spacing, bins=100, kde=True)
- plt.title("Latitude Spacing Distribution")
- plt.xlabel("Δlat (degrees)")
- plt.ylabel("Count")
- plt.tight_layout()
- plt.show()
-
- # Histogram: longitude spacing
- plt.figure(figsize=(10, 4))
- sns.histplot(lon_spacing, bins=100, kde=True)
- plt.title("Longitude Spacing Distribution")
- plt.xlabel("Δlon (degrees)")
- plt.ylabel("Count")
- plt.tight_layout()
- plt.show()
-
- # Density map (data space) — shift longitudes so the seam is at 180°
- lon_disp = lons_centered_for_display(df["lon"].values, center=180.0)
- plt.figure(figsize=(8, 6))
- plt.hexbin(lon_disp, df["lat"], gridsize=100, cmap="Reds", mincnt=1)
- plt.colorbar(label="Number of Points")
- plt.title("Grid Point Density (centered at 180°)")
- plt.xlabel("Longitude")
- plt.ylabel("Latitude")
- plt.tight_layout()
- plt.show()
-
-
-def plot_density_map(df: pd.DataFrame) -> None:
- """
- Plot a density map of grid points.
-
- Parameters
- ----------
- df : pd.DataFrame
- DataFrame containing lat/lon data.
- """
- lon_disp = lons_centered_for_display(df["lon"].values, center=180.0)
- plt.figure(figsize=(10, 6))
- plt.hexbin(lon_disp, df["lat"], gridsize=100, cmap="Reds", mincnt=1)
- plt.colorbar(label="Number of Points")
- plt.title("Grid Point Density (centered at 180°)")
- plt.xlabel("Longitude")
- plt.ylabel("Latitude")
- plt.tight_layout()
- plt.show()
-
-
-@app.command()
-def analyze(
- input_file: Annotated[
- Path,
- typer.Argument(
- exists=True,
- dir_okay=False,
- help="Input tomography file (.txt or .h5)"
- ),
- ],
-) -> None:
- """
- Analyze and plot tomography data from TXT or HDF5.
-
- Parameters
- ----------
- input_file : Path
- Input tomography file (.txt or .h5).
- """
- ext = input_file.suffix.lower()
- if ext == ".h5":
- df = load_hdf5_data(str(input_file))
- check_latlon_consistency(str(input_file))
- suffix = "(HDF5)"
- elif ext == ".txt":
- df = load_txt_ep_format(str(input_file))
- suffix = "(TXT)"
- else:
- raise typer.BadParameter("Unsupported file type. Use .txt or .h5")
-
- print(f"Loaded {len(df)} points.")
- print("Unique latitudes:", len(np.unique(df["lat"])))
- print("Unique longitudes:", len(np.unique(df["lon"])))
-
- lat_vals = np.sort(np.unique(df["lat"]))
- lon_vals = np.sort(np.unique(df["lon"]))
-
- print(f"Lat range: {df['lat'].min()} {df['lat'].max()}")
- print(f"Lon range: {df['lon'].min()} {df['lon'].max()}")
-
- check_grid_spacing(lat_vals, lon_vals)
-
- # Map view that handles dateline properly
- plot_spatial_distribution(df, suffix)
-
- # Spacing + data-space density
- plot_spacing_histograms(df)
-
- # If a depth column exists, show surface slice density
- if "depth" in df.columns:
- df_surface = df[df["depth"] == df["depth"].min()]
- if not df_surface.empty:
- plot_density_map(df_surface)
-
-
-if __name__ == "__main__":
- app()
diff --git a/tomography/tools/tomo_h5_downcast_to_f32.py b/tomography/tools/tomo_h5_downcast_to_f32.py
deleted file mode 100644
index 8e80b8a..0000000
--- a/tomography/tools/tomo_h5_downcast_to_f32.py
+++ /dev/null
@@ -1,133 +0,0 @@
-#!/usr/bin/env python3
-"""
-Downcast HDF5 tomography data from float64 to float32 with compression.
-
-This script converts vp, vs, and rho datasets from float64 to float32
-while recompressing the entire file with gzip compression and shuffle filter.
-"""
-
-from pathlib import Path
-from typing import Annotated
-
-import h5py
-import numpy as np
-import typer
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-TARGETS = {"vp", "vs", "rho"} # downcast only these
-
-
-def copy_attrs(src: h5py.Group | h5py.Dataset, dst: h5py.Group | h5py.Dataset) -> None:
- """
- Copy all attributes from source to destination object.
-
- Parameters
- ----------
- src : h5py.Group or h5py.Dataset
- Source HDF5 object.
- dst : h5py.Group or h5py.Dataset
- Destination HDF5 object.
- """
- for k, v in src.attrs.items():
- dst.attrs[k] = v
-
-
-def chunk_guess(shape: tuple) -> tuple | bool:
- """
- Guess appropriate chunk size for HDF5 dataset.
-
- Parameters
- ----------
- shape : tuple
- Shape of the dataset.
-
- Returns
- -------
- tuple or True
- Chunk shape for 2D datasets, True for automatic chunking for 1D.
- """
- # Good default for your 2D (1400x800) datasets; let h5py decide for 1D.
- if len(shape) == 2:
- return (256, 256)
- return True
-
-
-def downcast_hdf5(inp: str, outp: str, gzip_level: int = 4) -> None:
- """
- Convert HDF5 file with float64 to float32 downcast and recompression.
-
- Parameters
- ----------
- inp : str
- Input HDF5 file path.
- outp : str
- Output HDF5 file path.
- gzip_level : int, optional
- Gzip compression level (1-9). Default is 4.
- """
- with h5py.File(inp, "r") as fin, h5py.File(outp, "w") as fout:
- copy_attrs(fin, fout)
-
- def visit(name: str, obj: h5py.Group | h5py.Dataset) -> None:
- if isinstance(obj, h5py.Group):
- if name != "/":
- g = fout.require_group(name)
- copy_attrs(obj, g)
- elif isinstance(obj, h5py.Dataset):
- base = name.split("/")[-1]
- data = obj[()] # load to memory (fast path for recompress + retype)
- dtype = obj.dtype
- if (
- base in TARGETS
- and np.issubdtype(dtype, np.floating)
- and dtype != np.float32
- ):
- data = data.astype(np.float32, copy=False)
- dtype = np.float32
-
- ds = fout.create_dataset(
- name,
- data=data,
- compression="gzip",
- compression_opts=gzip_level,
- shuffle=True,
- chunks=chunk_guess(obj.shape),
- )
- copy_attrs(obj, ds)
-
- fin.visititems(lambda n, o: visit(n, o))
-
-
-@app.command()
-def main(
- input_file: Annotated[
- Path,
- typer.Argument(exists=True, dir_okay=False, help="Input HDF5 file path"),
- ],
- output_file: Annotated[
- Path,
- typer.Argument(dir_okay=False, help="Output HDF5 file path"),
- ],
- gzip: Annotated[
- int,
- typer.Option(min=1, max=9, help="Gzip compression level (1-9)"),
- ] = 4,
-) -> None:
- """
- Downcast vp/vs/rho to float32 and recompress with gzip+shuffle.
-
- Parameters
- ----------
- input_file : Path
- Input HDF5 file path.
- output_file : Path
- Output HDF5 file path.
- gzip : int, optional
- Gzip compression level (1-9). Default is 4.
- """
- downcast_hdf5(str(input_file), str(output_file), gzip_level=gzip)
-
-
-if __name__ == "__main__":
- app()
diff --git a/tomography/tools/tomo_in2h5.py b/tomography/tools/tomo_in2h5.py
deleted file mode 100644
index 01a4867..0000000
--- a/tomography/tools/tomo_in2h5.py
+++ /dev/null
@@ -1,275 +0,0 @@
-"""
-Convert ASCII tomography files to HDF5 format (with dtype/compression options).
-
-This script converts ASCII tomography files to HDF5 format. The input directory contains
-files with names like "surf_tomography_vp_elev0p25.in", "surf_tomography_vs_elev0p25.in" etc, where
-"vp", "vs", and "rho" are the velocity types and "0p25" is the elevation. The script reads the
-elevation values from the file names and ensures that they match across all velocity types. The
-output HDF5 file will contain groups for each elevation, with datasets for latitudes, longitudes,
-and the velocity types.
-
-Example usage:
-python tomo_in2h5.py /path/to/tomography_files output_name --out-dir /path/to/output_dir --dtype f32 --gzip 6
-
-This will convert the tomography files in /path/to/tomography_files to a file named output_name.h5
-in the same directory. If --out-dir is specified, the output file will be saved in that directory
-instead.
-"""
-
-import re
-from datetime import datetime
-from pathlib import Path
-from typing import Annotated
-
-import h5py
-import numpy as np
-import typer
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-
-def get_elevations_from_files(input_dir: Path) -> set[float]:
- """
- Extract unique elevation values from file names in input_dir.
-
- Parameters
- ----------
- input_dir : Path
- Path to directory containing tomography files.
-
- Returns
- -------
- set[float]
- Set of unique elevation values found in file names.
-
- Raises
- ------
- ValueError
- If no elevation files are found for any velocity type, or if elevations do not match
- across all velocity types.
-
- """
- vtypes = ["vp", "vs", "rho"]
- elev_pattern = re.compile(r"surf_tomography_(vp|vs|rho)_elev(-?\d+(?:p\d+)?)\.in")
- elevations_by_type: dict[str, set[float]] = {v: set() for v in vtypes}
- for filename in input_dir.glob("surf_tomography_*.in"):
- match = elev_pattern.match(filename.name)
- if match:
- vtype, elev_str = match.groups()
- elev = float(elev_str.replace("p", "."))
- elevations_by_type[vtype].add(elev)
- print(f"Found file: {filename.name}, vtype: {vtype}, elev: {elev}")
- if (
- not elevations_by_type["vp"]
- or not elevations_by_type["vs"]
- or not elevations_by_type["rho"]
- ):
- missing = [vtype for vtype, elevs in elevations_by_type.items() if not elevs]
- raise ValueError(
- f"No elevation files found for {', '.join(missing)} in {input_dir}"
- )
- if (
- elevations_by_type["vp"] != elevations_by_type["vs"]
- or elevations_by_type["vp"] != elevations_by_type["rho"]
- ):
- raise ValueError(
- "Elevations do not match across vp, vs, and rho:\n"
- f"vp: {sorted(elevations_by_type['vp'])}\n"
- f"vs: {sorted(elevations_by_type['vs'])}\n"
- f"rho: {sorted(elevations_by_type['rho'])}"
- )
- return elevations_by_type["vp"]
-
-
-def convert_ascii_to_hdf5(
- input_dir: Path,
- name: str,
- out_dir: Path | None = None,
- dtype_opt: str = "f64", # "f64" or "f32" for vp/vs/rho
- gzip_level: int = 4, # 1..9
- shuffle: bool = True,
-):
- """
- Convert ASCII tomography files to HDF5 format.
-
- Parameters
- ----------
- input_dir : Path
- Path to directory containing ASCII tomography files
- name : str
- Name for the output HDF5 file (without extension)
-
- out_dir : Path, optional
- Output directory for the HDF5 file (if unspecified, saves to input_dir/name.h5)
- dtype_opt : str, optional
- Data dtype for vp/vs/rho: "f64" (default) or "f32".
- gzip_level : int, optional
- Gzip compression level 1..9 (default 4).
- shuffle : bool, optional
- Enable shuffle filter (default is True).
-
- """
-
- vtypes = ["vp", "vs", "rho"]
- elevations = sorted(get_elevations_from_files(input_dir))
- if not elevations:
- raise ValueError(f"No valid tomography files found in {input_dir}")
- input_path = Path(input_dir)
-
- # Output path
- if out_dir is None:
- output_path = input_path / f"{name}.h5"
- else:
- Path(out_dir).mkdir(parents=True, exist_ok=True)
- output_path = Path(out_dir) / f"{name}.h5"
-
- # Dtype policy
- data_dtype = np.float64 if dtype_opt.lower() == "f64" else np.float32
- coord_dtype = np.float64 # keep coords precise
-
- with h5py.File(output_path, "w") as h5f:
- # file-level metadata to advertise decisions
- h5f.attrs["created"] = datetime.utcnow().isoformat() + "Z"
- h5f.attrs["generator"] = "tomo_in2h5.py"
- h5f.attrs["schema"] = "NZTomographyLevelStacked v1"
- h5f.attrs["data_dtype_vp_vs_rho"] = (
- "float32" if data_dtype == np.float32 else "float64"
- )
- h5f.attrs["coord_dtype_lat_lon"] = "float64"
- h5f.attrs["compression"] = f"gzip:{gzip_level}"
- h5f.attrs["shuffle"] = bool(shuffle)
-
- for elev in elevations:
- elev_file_str = (
- str(int(elev)) if elev == int(elev) else f"{elev:.2f}".replace(".", "p")
- )
- elev_group_name = str(int(elev)) if elev == int(elev) else f"{elev:.2f}"
- print(
- f"Processing elevation {elev_group_name} (files suffix: {elev_file_str})"
- )
- g = h5f.create_group(elev_group_name)
-
- # Load coords from rho file
- ref_file = input_path / f"surf_tomography_rho_elev{elev_file_str}.in"
- if not ref_file.exists():
- print(f"Warning: missing {ref_file}, skipping elevation {elev}")
- continue
-
- # Read header and coordinates
- with open(ref_file, "r") as f:
- nlat, nlon = map(int, f.readline().split())
- latitudes = np.array(
- [float(x) for x in f.readline().split()], dtype=coord_dtype
- )
- longitudes = np.array(
- [float(x) for x in f.readline().split()], dtype=coord_dtype
- )
-
- # coords - let h5py auto-chunk 1D arrays
- g.create_dataset(
- "latitudes",
- data=latitudes,
- compression="gzip",
- compression_opts=gzip_level,
- shuffle=shuffle,
- )
- g.create_dataset(
- "longitudes",
- data=longitudes,
- compression="gzip",
- compression_opts=gzip_level,
- shuffle=shuffle,
- )
-
- # fields
- for vtype in vtypes:
- filename = (
- input_path / f"surf_tomography_{vtype}_elev{elev_file_str}.in"
- )
- if not filename.exists():
- print(f"Warning: missing {filename}")
- continue
-
- # Use genfromtxt to skip header lines and load data directly
- arr = np.genfromtxt(filename, skip_header=3, dtype=data_dtype)
- arr = arr.reshape(nlat, nlon)
-
- dset = g.create_dataset(
- vtype,
- data=arr,
- compression="gzip",
- compression_opts=gzip_level,
- shuffle=shuffle,
- chunks=True,
- )
- dset.attrs["units"] = (
- "km/s" if vtype in ("vp", "vs") else "g/cm^3"
- ) # adjust if needed
- dset.attrs["dtype"] = (
- "float32" if data_dtype == np.float32 else "float64"
- )
-
- print(f"Wrote {vtype} at elevation {elev_group_name} to {output_path}")
- print(f"Done: {output_path}")
-
-
-@app.command()
-def convert_tomo_to_h5(
- input_dir: Annotated[
- Path,
- typer.Argument(
- exists=True,
- file_okay=False,
- help="Input directory containing ASCII tomography files",
- ),
- ],
- name: Annotated[
- str,
- typer.Argument(help="Base name for the output HDF5 file (without extension)"),
- ],
- out_dir: Annotated[
- Path | None,
- typer.Option(file_okay=False, help="Output directory for the HDF5 file"),
- ] = None,
- dtype: Annotated[
- str, typer.Option(help="Data dtype for vp/vs/rho: f32 (default) or f64")
- ] = "f32",
- gzip: Annotated[int, typer.Option(help="gzip level 1..9 (default 4)")] = 4,
- shuffle: Annotated[bool, typer.Option(help="Enable shuffle filter")] = True,
-):
- """
- Convert ASCII tomography files to HDF5 format.
-
- This tool consolidates multiple ASCII tomography files into a single,
- structured HDF5 file. It handles files with various elevation values
- and velocity parameters (vp, vs, rho) and ensures they are consistent.
-
- Parameters
- ----------
- input_dir : Path
- Path to directory containing ASCII tomography files.
- name : str
- Name for the output HDF5 file (without extension).
- out_dir : Path, optional
- Output directory for the HDF5 file (if unspecified, saves to input_dir/name.h5).
- dtype : str, optional
- Data dtype for vp/vs/rho: "f64" (default) or "f32".
- gzip : int, optional
- Gzip compression level 1..9 (default 4).
- shuffle : bool, optional
- Enable shuffle filter (default is True).
-
- """
-
- convert_ascii_to_hdf5(
- input_dir,
- name,
- out_dir,
- dtype_opt=dtype,
- gzip_level=gzip,
- shuffle=shuffle,
- )
-
-
-if __name__ == "__main__":
- app()
diff --git a/tomography/tools/visualisation/tomo_3dviewer.py b/tomography/tools/visualisation/tomo_3dviewer.py
deleted file mode 100644
index dbe6319..0000000
--- a/tomography/tools/visualisation/tomo_3dviewer.py
+++ /dev/null
@@ -1,561 +0,0 @@
-"""
-3D Tomography Viewer for HDF5-based velocity models.
-
-This script provides an interactive 3D visualization of tomography slices
-stored in an HDF5 file. It uses PyVista and PyQt for the user interface
-and rendering. Users can toggle the visibility of different elevation layers,
-overlay original data points from a text file, and inspect the model from
-various angles.
-
-The script is run from the command line and accepts various options to
-customize the visualization, such as selecting the scalar field, defining
-latitude/longitude ranges, and setting color map limits.
-"""
-
-import sys
-from pathlib import Path
-from typing import Annotated, Optional
-
-import h5py
-import numpy as np
-import pandas as pd
-import pyvista as pv
-import typer
-from pyvista.plotting.camera import Camera
-from pyvistaqt import QtInteractor
-from qtpy.QtCore import Qt
-from qtpy.QtWidgets import (
- QApplication,
- QCheckBox,
- QDockWidget,
- QLabel,
- QMainWindow,
- QScrollArea,
- QVBoxLayout,
- QWidget,
-)
-
-from qcore import cli
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-DEFAULT_CAMERA_ANGLE = (-0.018774705983602088, 0.8058896333404005, 0.5917680367252902)
-
-
-def is_number(s: str) -> bool:
- """Check if a string can be converted to a float.
-
- Parameters
- ----------
- s : str
- The input string.
-
- Returns
- -------
- bool
- True if the string can be converted to a float, False otherwise.
- """
- try:
- float(s)
- return True
- except ValueError:
- return False
-
-
-def read_ep_txt(txt_path: Path) -> pd.DataFrame:
- """Read an EP-style tomography text file into a pandas DataFrame.
-
- The file is expected to have a specific format with columns for velocity
- parameters and coordinates. This function handles file parsing, column
- naming, and coordinate adjustments (elevation sign and longitude range).
-
- Parameters
- ----------
- txt_path : Path
- Path to the input text file.
-
- Returns
- -------
- pd.DataFrame
- A DataFrame containing the tomography data with standardized column names.
- """
- col_names = [
- "vp",
- "vp_o_vs",
- "vs",
- "rho",
- "sf_vp",
- "sf_vp_o_vs",
- "x",
- "y",
- "elevation",
- "lat",
- "lon",
- ]
- df = pd.read_csv(txt_path, sep=r"\s+", skiprows=2, names=col_names, engine="python")
- df["elevation"] = -1 * df["elevation"]
- df["lon"] = np.where(df["lon"] < 0, df["lon"] + 360, df["lon"])
- return df
-
-
-def load_stacked_slices(
- h5_path: Path,
- scalar: str = "vp",
- lat_range: Optional[tuple[float, float]] = None,
- lon_range: Optional[tuple[float, float]] = None,
-) -> tuple[list[str], np.ndarray, np.ndarray, np.ndarray]:
- """Load and stack tomography slices from an HDF5 file.
-
- Reads data for specified elevation layers from an HDF5 file, optionally
- cropping it to a given latitude and longitude range.
-
- Parameters
- ----------
- h5_path : Path
- Path to the HDF5 file.
- scalar : str, optional
- The scalar field to load (e.g., "vp", "vs"), by default "vp".
- lat_range : tuple of float, optional
- The latitude range (min, max) to crop the data. If None, all latitudes
- are used. By default None.
- lon_range : tuple of float, optional
- The longitude range (min, max) to crop the data. If None, all latitudes
- are used. By default None.
-
- Returns
- -------
- tuple
- A tuple containing:
- - keys (list[str]): Sorted list of elevation keys.
- - lon (np.ndarray): Array of longitudes within the specified range.
- - lat (np.ndarray): Array of latitudes within the specified range.
- - scalar_values np.ndarray: A 3d array containing scalar data of shape (n_elevations, n_latitudes, n_longitudes).
- """
- with h5py.File(h5_path, "r") as f:
- keys = sorted((k for k in f.keys() if is_number(k)), key=float, reverse=True)
- sample = f[keys[0]]
- lat_full = sample["latitudes"][:]
- lon_full = sample["longitudes"][:]
-
- lat_mask = (
- (lat_full >= lat_range[0]) & (lat_full <= lat_range[1])
- if lat_range
- else np.ones_like(lat_full, bool)
- )
- lon_mask = (
- (lon_full >= lon_range[0]) & (lon_full <= lon_range[1])
- if lon_range
- else np.ones_like(lon_full, bool)
- )
-
- lat = lat_full[lat_mask]
- lon = lon_full[lon_mask]
-
- scalar_values = []
- for k in keys:
- data = f[k][scalar][:][np.ix_(lat_mask, lon_mask)]
- scalar_values.append(data)
- return keys, lon, lat, np.array(scalar_values)
-
-
-def make_flat_surfaces(
- elevations: list[str],
- lon: np.ndarray,
- lat: np.ndarray,
- scalar_values: np.ndarray,
- gap: float = 0.1,
-) -> list[tuple[pv.StructuredGrid, float, float, float]]:
- """Create a list of PyVista surfaces from stacked 2D data.
-
- Each slice in the data cube is converted into a flat surface mesh,
- vertically stacked with a specified gap.
-
- Parameters
- ----------
- elevations : list[str]
- List of elevation values for each slice.
- lon : np.ndarray
- 1D array of longitudes.
- lat : np.ndarray
- 1D array of latitudes.
- scalar_values : np.ndarray
- 3D data array with shape (elevation, lat, lon).
- gap : float, optional
- Vertical gap between the stacked surfaces, by default 0.1.
-
-
- Returns
- -------
- tuple
- A tuple containing:
- - grids (dict): A dictionary mapping elevations to their corresponding PyVista surface meshes.
- Each key is a float elevation value, and the value is a PyVista StructuredGrid object.
- - global_min (float): The minimum value in the entire data cube, used for consistent color mapping.
- - global_max (float): The maximum value in the entire data cube, used for consistent color mapping.
-
- """
-
- grid_dict = {}
- global_min = np.min(scalar_values)
- global_max = np.max(scalar_values)
-
- for iz in range(len(elevations)):
- z = -gap * iz
- xs, ys = np.meshgrid(lon, lat, indexing="xy")
- zs = np.full_like(xs, z)
- surf = pv.StructuredGrid(xs, ys, zs).extract_surface()
- surf["values"] = scalar_values[iz].ravel(order="C")
- grid_dict[float(elevations[iz])] = surf
-
- return grid_dict, global_min, global_max
-
-
-class TomoApp(QMainWindow):
- """Main application window for the 3D tomography viewer.
-
- This class encapsulates the Qt-based GUI, including the PyVista plotter
- and layer visibility controls.
-
-
- Parameters
- ----------
- title : str
- Title for the window.
- scalar_name : str
- Name of the scalar field being displayed.
- grid_dict : dict[pv.StructuredGrid]
- Dictionary mapping elevations to PyVista StructuredGrid objects.
- clim : tuple, optional
- Color map limits (min, max) to override the defaults. By default None.
- points_by_elevation : dict, optional
- Dictionary mapping elevations to point data to be overlaid. By default None.
- debug : bool, optional
- Flag to enable debug printing. By default False.
-
- """
-
- def __init__(
- self,
- title: str,
- scalar_name: str,
- grid_dict: dict[pv.StructuredGrid],
- clim: tuple[float, float],
- points_by_elevation: Optional[dict[float, np.ndarray]] = None,
- debug: bool = False,
- ):
- """Initialize the Tomography 3D Viewer application."""
-
- super().__init__()
- self.debug = debug
-
- # Find the topmost grid (highest elevation)
- top_grid = grid_dict[max(grid_dict.keys())]
- # Dynamically determine focal point from the center of the top grid
- self.focal_point = top_grid.center
-
- # Store camera settings relative to the focal point
- self.camera_position = (
- self.focal_point[0] - 1.63,
- self.focal_point[1] - 46.57,
- self.focal_point[2] + 63.37,
- )
- self.view_up = DEFAULT_CAMERA_ANGLE
-
- self.plotter_widget = QtInteractor(self)
-
- self.setWindowTitle(f"Tomography 3D Viewer - {title}")
-
- title_label = QLabel(f"{title}- {scalar_name}")
- title_label.setAlignment(Qt.AlignCenter)
- title_label.setStyleSheet("font-weight: bold; font-size: 14px;")
-
- main_layout = QVBoxLayout()
- main_layout.addWidget(title_label)
- main_layout.addWidget(self.plotter_widget)
- central_widget = QWidget()
- central_widget.setLayout(main_layout)
- self.setCentralWidget(central_widget)
- self.actors = []
- self.slice_actor = None
- self.point_actors = []
- self.points_by_elevation = points_by_elevation
-
- self.toggle_panel = QWidget()
- layout = QVBoxLayout(self.toggle_panel)
- # layout.addWidget(QLabel(self.windowTitle()))
- layout.addWidget(QLabel("Layer Visibility"))
-
- scroll = QScrollArea()
- scroll.setWidgetResizable(True)
- scroll.setWidget(self.toggle_panel)
-
- dock = QDockWidget("Layers", self)
- dock.setWidget(scroll)
- dock.setAllowedAreas(Qt.LeftDockWidgetArea | Qt.RightDockWidgetArea)
- self.addDockWidget(Qt.LeftDockWidgetArea, dock)
-
- for i, elevation in enumerate(
- sorted(grid_dict.keys())[::-1]
- ): # reverse order for top-down view
- grid = grid_dict[elevation]
- clim = clim
- actor = self.plotter_widget.add_mesh(
- grid,
- scalars="values",
- cmap="hot_r",
- opacity=1.0,
- clim=clim,
- show_scalar_bar=False,
- name=f"{elevation}km",
- lighting=False,
- )
- self.actors.append(actor)
-
- cb = QCheckBox(f"{elevation:.0f} km")
- cb.setChecked(True)
-
- # Add optional point overlay for this elevation if provided
- if self.points_by_elevation and elevation in self.points_by_elevation:
- z = -i * 0.1 + 0.01 # slightly above the surface
- pts = self.points_by_elevation[elevation]
- pdata = pv.PolyData(
- np.column_stack((pts[:, 1], pts[:, 0], np.full(len(pts), z)))
- )
- pdata["values"] = pts[:, 2] if pts.shape[1] > 2 else np.zeros(len(pts))
- if self.debug:
- print(pdata["values"])
-
- p_actor = self.plotter_widget.add_mesh(
- pdata,
- scalars="values",
- cmap="hot_r",
- point_size=4,
- render_points_as_spheres=True,
- show_scalar_bar=False,
- clim=clim,
- )
- self.point_actors.append((elevation, p_actor))
- if elevation == max(self.points_by_elevation):
- p_actor.SetVisibility(True)
- else:
- p_actor.SetVisibility(False)
-
- cb.stateChanged.connect(
- lambda state, a=actor: self.set_visibility(a, state)
- )
- layout.addWidget(cb)
-
- self.plotter_widget.add_scalar_bar(title=scalar_name.upper(), n_labels=5)
- self.plotter_widget.add_axes(
- xlabel="Longitude (°E)",
- ylabel="Latitude (°N)",
- zlabel="Elevation",
- label_size=(0.4, 0.4),
- )
- self.plotter_widget.camera_position = "iso"
- self.plotter_widget.set_scale(xscale=1, yscale=1, zscale=0.05)
- self.plotter_widget.renderer.interpolate_before_map = True
-
- # 🧭 Set custom camera view
- self.reset_camera()
-
- self.plotter_widget.add_key_event("v", self.toggle_all)
- self.plotter_widget.add_key_event("s", self.toggle_slice)
- self.plotter_widget.add_key_event("r", self.reset_camera)
-
- self.plotter_widget.renderer.GetActiveCamera().AddObserver(
- "ModifiedEvent", self.on_camera_modified
- )
-
- def set_topmost_point_visibility(self) -> None:
- """Control visibility of point overlays.
-
- Ensures that only the points corresponding to the topmost visible
- elevation layer are displayed.
- """
- if not self.points_by_elevation:
- return
- checkboxes = self.toggle_panel.findChildren(QCheckBox)
- visible_elevations = [
- float(cb.text().split()[0]) for cb in checkboxes if cb.isChecked()
- ]
- if self.debug:
- print(visible_elevations)
- if visible_elevations:
- top = max(visible_elevations) #
- for d, a in self.point_actors:
- is_visible = d == top
- a.SetVisibility(is_visible)
- if is_visible:
- if self.debug:
- print(f"👁️ Top visible elevation: {top} km")
- pdata = a.GetMapper().GetInputAsDataSet()
- if "values" in pdata.point_data:
- if self.debug:
- print("📊 Scalar values:", pdata["values"])
- else:
- if self.debug:
- print("⚠️ 'values' not found in point data")
-
- def set_visibility(self, actor: pv.Actor, state: int) -> None:
- """Set the visibility of a layer and update point overlays.
-
- This is a slot connected to the layer visibility checkboxes.
-
- Parameters
- ----------
- actor : pv.Actor
- The main surface actor for the layer.
-
- state : int
- The state of the checkbox (Qt.Checked or Qt.Unchecked).
-
- """
- actor.SetVisibility(state == Qt.Checked)
- self.set_topmost_point_visibility()
- self.plotter_widget.render()
-
- def toggle_all(self) -> None:
- """Toggle the visibility of all layers."""
- for actor in self.actors:
- actor.SetVisibility(not actor.GetVisibility())
- self.plotter_widget.render()
-
- def toggle_slice(self) -> None:
- """Toggle a cross-section slice view."""
- if self.slice_actor:
- self.plotter_widget.remove_actor(self.slice_actor)
- self.slice_actor = None
- else:
- grid = self.actors[len(self.actors) // 2].GetMapper().GetInputAsDataSet()
- bounds = grid.bounds
- x_center = 0.5 * (bounds[0] + bounds[1])
- sliced = grid.slice(normal="x", origin=(x_center, 0, 0))
- self.slice_actor = self.plotter_widget.add_mesh(
- sliced, color="white", line_width=2
- )
- self.plotter_widget.render()
-
- def reset_camera(self) -> None:
- """Reset the camera to a predefined position and orientation."""
- camera = self.plotter_widget.renderer.GetActiveCamera()
- camera.SetPosition(self.camera_position)
- camera.SetFocalPoint(self.focal_point)
- camera.SetViewUp(self.view_up)
- self.plotter_widget.renderer.ResetCameraClippingRange()
- self.plotter_widget.render()
-
- def on_camera_modified(self, caller: Camera, event: str) -> None:
- """Print camera parameters when modified (for debugging).
-
- Parameters
- ----------
- caller : Camera
- The camera object that triggered the event.
- event : str
- The event name (e.g., "ModifiedEvent").
- """
- if not self.debug:
- return
-
- cam = self.plotter_widget.renderer.GetActiveCamera()
- print("📸 Camera changed:")
- print(" Position :", cam.GetPosition())
- print(" FocalPoint :", cam.GetFocalPoint())
- print(" ViewUp :", cam.GetViewUp())
-
-
-@cli.from_docstring(app)
-def launch_viewer(
- h5file: Annotated[Path, typer.Argument(help="HDF5 file with tomography slices")],
- scalar: Annotated[
- str, typer.Option(help="Scalar field to show", show_default=True)
- ] = "vp",
- gap: Annotated[float, typer.Option(help="Vertical gap between slices")] = 0.2,
- lat_range: Annotated[
- Optional[tuple[float, float]], typer.Option(help="Latitude range [min max]")
- ] = None,
- lon_range: Annotated[
- Optional[tuple[float, float]], typer.Option(help="Longitude range [min max]")
- ] = None,
- clim: Annotated[
- Optional[tuple[float, float]], typer.Option(help="Color range [min max]")
- ] = None,
- txt: Annotated[
- Optional[Path],
- typer.Option(help="Optional EP-style TXT input file for original points"),
- ] = None,
- debug: Annotated[bool, typer.Option(help="Enable debug print statements")] = False,
-) -> None:
- """Launch the interactive tomography viewer with optional overlays and controls.
-
- This function serves as the main entry point for the command-line application.
- It parses arguments, loads data, and initializes the Qt application and
- the TomoApp viewer window.
-
- Parameters
- ----------
- h5file : Path
- Path to the HDF5 file with tomography slices.
- scalar : str, optional
- Scalar field to show (e.g., "vp", "vs"), by default "vp".
- gap : float, optional
- Vertical gap between slices, by default 0.2.
- lat_range : tuple of float, optional
- Latitude range (min, max) to display, by default None.
- lon_range : tuple of float, optional
- Longitude range (min, max) to display, by default None.
- clim : tuple of float, optional
- Color range (min, max) for the scalar bar, by default None.
- txt : Path, optional
- Optional EP-style TXT input file for original points overlay, by default None.
- debug : bool, optional
- Enable debug print statements, by default False.
- """
- elevations, lon, lat, scalar_values = load_stacked_slices(
- h5file, scalar=scalar, lat_range=lat_range, lon_range=lon_range
- )
- grid_dict, vmin, vmax = make_flat_surfaces(
- elevations, lon, lat, scalar_values, gap=gap
- )
-
- points_by_elevation = None
- if txt:
- df = read_ep_txt(txt)
- df_filtered = df[np.any(np.isclose(df['elevation'].values[:, None], [float(d) for d in elevations], atol=1e-3), axis=1)]
- if debug:
- print(df_filtered.head())
- points_by_elevation = {
- float(d): df_filtered[
- np.isclose(df_filtered["elevation"], float(d), atol=1e-3)
- ][["lat", "lon", scalar]].values
- for d in elevations
- }
-
- app_qt = QApplication.instance() or QApplication(sys.argv)
-
- # validate clim input
- if clim is None:
- clim = (vmin, vmax)
- else:
- if len(clim) != 2 or not all(isinstance(c, (int, float)) for c in clim):
- raise ValueError("clim must be a tuple of two numeric values (min, max).")
- clim = tuple(clim)
-
- window = TomoApp(
- h5file.stem,
- scalar,
- grid_dict,
- clim,
- points_by_elevation=points_by_elevation,
- debug=debug,
- )
- window.show()
- # Ensure clean shutdown
- app_qt.lastWindowClosed.connect(app_qt.quit)
-
- sys.exit(app_qt.exec())
-
-
-if __name__ == "__main__":
- app()
diff --git a/tomography/tools/visualisation/tomo_boundary.py b/tomography/tools/visualisation/tomo_boundary.py
deleted file mode 100644
index 909b3a5..0000000
--- a/tomography/tools/visualisation/tomo_boundary.py
+++ /dev/null
@@ -1,100 +0,0 @@
-"""
-This script extracts the geographical boundary from a tomography HDF5 file
-and saves it as a GeoJSON file. It computes the convex hull of the grid points
-defined by the latitude and longitude arrays in the HDF5 file to determine
-the boundary polygon.
-
-The script is intended to be run from the command line, taking the path to an
-HDF5 file as input.
-"""
-
-from pathlib import Path
-from typing import Annotated
-
-import geojson
-import h5py
-import numpy as np
-import typer
-from shapely.geometry import MultiPoint
-
-from qcore import cli
-
-app = typer.Typer(pretty_exceptions_enable=False)
-
-
-@cli.from_docstring(app)
-def extract_single_boundary_geojson(
- hdf5_file_path: Annotated[
- Path,
- typer.Argument(
- help="Path to the HDF5 file (e.g., 2020_NZ.h5)",
- exists=True,
- file_okay=True,
- dir_okay=False,
- resolve_path=True,
- ),
- ],
-):
- """Extracts the convex hull boundary from an HDF5 file and saves it as GeoJSON.
-
- This function reads latitude and longitude grids from the first group in the
- HDF5 file, calculates the convex hull of the grid's edge points, and
- writes the resulting polygon to a GeoJSON file. The output file will have
- the same name as the input file but with a .geojson extension.
-
- Parameters
- ----------
- hdf5_file_path : Path
- The path to the input HDF5 file containing latitude and longitude grids.
-
- Raises
- ------
- RuntimeError
- If the convex hull operation does not result in a Polygon.
- """
- with h5py.File(hdf5_file_path, "r") as f:
- # Use the first available group to extract lat/lon
- first_group_name = next(iter(f.keys()))
- group = None
- for key in f.keys():
- if 'latitudes' in f[key] and 'longitudes' in f[key]:
- group = f[key]
- break
- if group is None:
- raise RuntimeError("No group with 'latitudes' and 'longitudes' found in HDF5 file.")
-
- lats = group["latitudes"][:]
- lons = group["longitudes"][:]
-
- # Construct edge points only
- lon_grid, lat_grid = np.meshgrid(lons, lats)
- edge_coords = []
-
- # Top and bottom rows
- edge_coords += list(zip(lon_grid[0, :], lat_grid[0, :])) # top
- edge_coords += list(zip(lon_grid[-1, :], lat_grid[-1, :])) # bottom
-
- # Left and right columns (excluding corners)
- edge_coords += list(zip(lon_grid[1:-1, 0], lat_grid[1:-1, 0])) # left
- edge_coords += list(zip(lon_grid[1:-1, -1], lat_grid[1:-1, -1])) # right
-
- # Compute convex hull
- hull = MultiPoint(edge_coords).convex_hull
-
- if hull.geom_type != "Polygon":
- raise RuntimeError("Convex hull did not produce a polygon.")
-
- coords = [list(p) for p in hull.exterior.coords]
- polygon = geojson.Polygon([coords])
- feature = geojson.Feature(geometry=polygon, properties={})
- feature_collection = geojson.FeatureCollection([feature])
-
- output_path = hdf5_file_path.with_suffix(".geojson")
- with open(output_path, "w") as geojson_file:
- geojson.dump(feature_collection, geojson_file, indent=2)
-
- print(f"Boundary GeoJSON saved to: {output_path}")
-
-
-if __name__ == "__main__":
- app()