Release Volume Sampler: Preparational Steps¶

This notebook lets you run each preparational step interactively.

In [1]:
# Imports
import os
import numpy as np
import shutil
import rasterio

# Change to the src directory
# This is necessary to ensure that the script runs in the correct context
# and can find the necessary modules and files.
os.chdir("/home/ebr/projects/release-volume-sampler/src")

The history saving thread hit an unexpected error (OperationalError('database or disk is full')).History will not be written to the database.
In [2]:
from rvsampler.preprocess import truncate_positive_values, slope, aspect
from rvsampler.slope_analysis import SlopeAnalysis
from rvsampler.triangulate import Triangulate, Triangulation
from rvsampler.release_volume_sampler import RecursiveReleaseAnalysis
from rvsampler.database_handler import VolumeDatabaseHandler
from rvsampler.cluster import ClusterAnalysis
from rvsampler.cumprobs_by_triangle import caclulate_cumulative_probabilities
from rvsampler.utils import create_dir
from rvsampler.set_logg import setup_logger

2025-08-07 07:35:22.718683: I tensorflow/core/util/port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2025-08-07 07:35:22.720253: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:32] Could not find cuda drivers on your machine, GPU will not be used.
2025-08-07 07:35:22.724627: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:32] Could not find cuda drivers on your machine, GPU will not be used.
2025-08-07 07:35:22.736266: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:477] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
E0000 00:00:1754552122.755593 1522505 cuda_dnn.cc:8310] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
E0000 00:00:1754552122.761327 1522505 cuda_blas.cc:1418] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
2025-08-07 07:35:22.782590: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F AVX512_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
In [3]:
# Set region and directories
region = "messina_20250806"  # Change as needed
rootdir = "/home/ebr/projects/release-volume-sampler"  # Change as needed
rundir = os.path.join(rootdir, 'generated', region)
logger = setup_logger("preparational", rundir)
create_dir(rundir, logger)
logger.info(f"Running preparational notebook for region {region} in {rundir}")

2025-08-07 07:35:26,101 - preparational - INFO - Running preparational notebook for region messina_20250806 in /home/ebr/projects/release-volume-sampler/generated/messina_20250806
In [ ]:
# Configuration
config = {
    "rundir": rundir,
    "singularity_image": os.path.join(rootdir,"images", "grass.sif"),
    "bathyfile": os.path.join(rootdir,'input', "bathy", "messina_001/bathy_truncated.tif"),
    "soilregions_filename": os.path.join(rootdir,'input', 'soilparams','regions.tif'),
    "soil_parameters_filename": os.path.join(rootdir,'input', 'soilparams','params.json'),
    "logger": logger,
    "slopeunitfile": None, #os.path.join(rootdir, 'input', 'slopeunits', 'slumap.tif'),
}
logger.info(f"Configuration: {config}")

# Initialization
def initialize(rundir, bathyfile, soilregions_filename, soil_parameters_filename, singularity_image, logger, slopeunitfile=None):
    logfile = os.path.join(rundir, "preparational_external_software.txt")
    shutil.copy(soilregions_filename, os.path.join(rundir, "soilregions.tif"))
    shutil.copy(soil_parameters_filename, os.path.join(rundir, "soilparams.json"))
    if slopeunitfile is not None:
        shutil.copy(slopeunitfile, os.path.join(rundir, "slumap.tif"))
    logger.info("Verifying that input bathymetri is logitude-latitude.")
    with rasterio.open(bathyfile) as src:
        assert src.crs.is_geographic, "The raster is not in a longitude-latitude coordinate system (geographic CRS)."
    logger.info(f"Copy {bathyfile} to {rundir}.")
    outfile = "bathy"    
    shutil.copy(bathyfile, os.path.join(rundir, f"{outfile}.tif"))
    outfile = truncate_positive_values(rundir, singularity_image, outfile, logfile) 
    slope(bathyfile, output_dir=rundir, logfile=logfile)
    aspect(bathyfile, output_dir=rundir, logfile=logfile)
    return rundir

# Run initialization
initialize(**config)

Execute Slope Analysis¶

In [ ]:
def execute_slope_analysis(rundir):
    sa = SlopeAnalysis(rundir, slopefile="slope.tif")
    quantiles = [0.01, 0.1, 0.5, 0.9, 0.99]
    sa.compute_quantiles(quantiles, write_fos=True, write_ky=True)
    fos_thresholds = np.linspace(0, 2, num=50)
    sa.compute_cumulative(fos_thresholds, feature_name="logfos", write=True)
    ky_thresholds = np.linspace(-3,1, num=50)
    sa.compute_cumulative(ky_thresholds, feature_name="logky", write=True)

execute_slope_analysis(rundir)

Triangulate domain¶

In [ ]:
config = {
    "rundir": rundir,
    "bathyfile": "bathy_truncated.tif",
    "utm_epsg_code": 32633,
    "resolution": (110, 110),
    "slopeunitfile": os.path.join(rundir, "slumap.tif"),
}
optimization_params = {
    "num_iterations": 400,
    "batch_size": 3000,
    "shape_weight": 5e1,
    "area_weight": 5e-11,
    "elevation_weight": 1e-2
}
triang = Triangulate(**config)
triang.fit(**optimization_params)
triang.write_to_file()
triang.plot_triangulation(output_file="triangulation.png")
In [ ]:
# Reload triangulation and compute aditional quantities.
triang = Triangulation(rundir)

triang.initialize_triangle_properties()
triang.poly_slopes(filename="poly_slopes.npy")
triang.slopeunits_to_triangles(filename="slopeunits_to_triangles.npy")
cumulative_dir = os.path.join(rundir, "slope_analysis", "fos", "cumulative")
outfile_name = "cumulative_fos.npz"
triang.create_lookuptable(cumulative_dir, outfile_name)

Sample Release Volumes¶

In [ ]:
with VolumeDatabaseHandler(rundir) as db:
    db.initialize_db()

config = {
    "rundir": rundir,
    "mesh_path": os.path.join(rundir, "triangulation", "triangulation.vtk"),
    "cumprob_logfos_path": os.path.join(rundir, "slope_analysis", "fos", "cumulative", "cumulative_fos.npz"),
    "utm_epsg_code": 32633,
}

run_config = {
    "fos_threshold": 1.6,
    "recursive_probability_threshold": 1.e-4,
    "seed_triangle_probability_threshold": 0.05,
    "max_workers":10,
    "max_n_seed_triangles": 2000,
    "use_slopeunits": False,
    "max_n_slopeunits": 10000,
    "max_n_simultaneous": 2,
    "pairs_distance_threshold": 1000,
}
analysis = RecursiveReleaseAnalysis(**config)
analysis.run(**run_config)

with VolumeDatabaseHandler(rundir) as db:
    assert db.test_no_duplicate_triangles_in_released()
In [ ]:
# Assign volume, thickness, tsunami potential ratio and no2d to all volumes
volume_factor = 2.0  # Factor to scale the thickness of the release volumes. 

with VolumeDatabaseHandler(rundir) as db:
    db.assign_volume_features(volume_factor=volume_factor)
In [7]:
with VolumeDatabaseHandler(rundir) as volumes_db:
    volumes_db.plot_distribution(seed_prob="p_fos_seed")

2025-08-07 07:37:20,044 - db_handler - INFO - Database connection initialized.
2025-08-07 07:37:44,737 - db_handler - INFO - Database connection closed.
No description has been provided for this image
In [8]:
with VolumeDatabaseHandler(rundir) as volumes_db:
    volumes_db.plot_release_density_plots(seed_prob="p_fos_seed")

2025-08-07 07:37:45,716 - db_handler - INFO - Database connection initialized.
2025-08-07 07:37:45,718 - db_handler - INFO -  Create spatial release density figure.
2025-08-07 07:38:53,207 - db_handler - INFO - Aggregated triangles: 0
2025-08-07 07:38:57,686 - db_handler - INFO - Aggregated triangles: 500
2025-08-07 07:39:01,954 - db_handler - INFO - Aggregated triangles: 1000
2025-08-07 07:39:06,237 - db_handler - INFO - Aggregated triangles: 1500
2025-08-07 07:39:10,586 - db_handler - INFO - Aggregated triangles: 2000
2025-08-07 07:39:14,839 - db_handler - INFO - Aggregated triangles: 2500
2025-08-07 07:39:19,189 - db_handler - INFO - Aggregated triangles: 3000
2025-08-07 07:39:23,492 - db_handler - INFO - Aggregated triangles: 3500
2025-08-07 07:39:27,850 - db_handler - INFO - Aggregated triangles: 4000
2025-08-07 07:39:32,152 - db_handler - INFO - Aggregated triangles: 4500
2025-08-07 07:39:36,464 - db_handler - INFO - Aggregated triangles: 5000
2025-08-07 07:39:40,798 - db_handler - INFO - Aggregated triangles: 5500
2025-08-07 07:39:45,786 - db_handler - INFO - Aggregated triangles: 6000
2025-08-07 07:39:50,106 - db_handler - INFO - Aggregated triangles: 6500
2025-08-07 07:39:54,407 - db_handler - INFO - Aggregated triangles: 7000
2025-08-07 07:39:58,697 - db_handler - INFO - Aggregated triangles: 7500
2025-08-07 07:40:02,999 - db_handler - INFO - Aggregated triangles: 8000
2025-08-07 07:40:07,400 - db_handler - INFO - Aggregated triangles: 8500
2025-08-07 07:40:11,692 - db_handler - INFO - Aggregated triangles: 9000
2025-08-07 07:40:16,020 - db_handler - INFO - Aggregated triangles: 9500
2025-08-07 07:40:20,423 - db_handler - INFO - Aggregated triangles: 10000
2025-08-07 07:40:25,105 - db_handler - INFO - Aggregated triangles: 10500
2025-08-07 07:40:29,698 - db_handler - INFO - Aggregated triangles: 11000
2025-08-07 07:40:34,111 - db_handler - INFO - Aggregated triangles: 11500
2025-08-07 07:40:38,469 - db_handler - INFO - Aggregated triangles: 12000
2025-08-07 07:40:42,805 - db_handler - INFO - Aggregated triangles: 12500
2025-08-07 07:40:47,132 - db_handler - INFO - Aggregated triangles: 13000
2025-08-07 07:40:51,705 - db_handler - INFO - Aggregated triangles: 13500
2025-08-07 07:40:56,236 - db_handler - INFO - Aggregated triangles: 14000
2025-08-07 07:41:00,786 - db_handler - INFO - Database connection closed.
No description has been provided for this image

Cluster release volumes¶

In [ ]:
def cluster_release_volumes(rundir):
    config = {
        "rundir": rundir,
        "n_clusters": 100,
        "random_state": 0,
        "batch_size": 1000,
        "feature_columns": ['area', 'no2d', 'mean_elevation', 'mean_northing', 'mean_easting'],
        "columns_to_scale": ['area', 'no2d', 'mean_elevation', 'mean_northing', 'mean_easting'],
        "weights": {
            'area': 1.0,
            'no2d': 1.0,
            'mean_elevation': 1.0,
            'mean_northing': 10.0,
            'mean_easting': 10.0
        },
    }
    cluster_analysis = ClusterAnalysis(**config)
    cluster_analysis.fit()
    cluster_analysis.write_to_database()
    cluster_analysis.find_representatives()
    cluster_analysis.close()

cluster_release_volumes(rundir)
In [ ]:
# Write cluster representatives to CSV and Rasters
def write_volumes(rundir):
    filter_config = {
        "tsunami_potential_ratio_threshold": 0.,
        "max_rasters": 100,
        "raster_driver": 'AAIGrid', #AAIGrid or GTiff
    }
    with VolumeDatabaseHandler(rundir) as volumes_db:
        volumes_db.write_volumes_to_csv(max_rasters=filter_config['max_rasters'])
        volumes_db.write_volumes_to_rasters(**filter_config)

write_volumes(rundir)
In [4]:
import numpy as np
import matplotlib.pyplot as plt
import rasterio

rasters_dir = os.path.join(rundir, "volumes/rasters")
bathy_path = os.path.join(rundir,  "bathy.tif")

# Read bathymetry
with rasterio.open(bathy_path) as src:
    bathy = src.read(1)
    bounds = src.bounds
    extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]

plt.figure(figsize=(10, 8))
plt.imshow(bathy, cmap='Blues', extent=extent, origin='upper')
plt.title("Release Volumes")
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# Plot all .tif rasters as contours
for fname in os.listdir(rasters_dir):
    if fname.endswith(".asc") or fname.endswith(".tif"):
        fpath = os.path.join(rasters_dir, fname)
        with rasterio.open(fpath) as src:
            data = src.read(1)
            bounds = src.bounds
            extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
            # Mask no-data values if needed
            data = np.ma.masked_equal(data, src.nodata)
            # Generate coordinates
            x = np.linspace(bounds.left, bounds.right, data.shape[1])
            y = np.linspace(bounds.bottom, bounds.top, data.shape[0])
            X, Y = np.meshgrid(x, y)
            # Plot contour for nonzero values
            if np.any(data > 0):
                plt.contour(X, Y, data, levels=[0], colors='r', linewidths=1, alpha=0.5)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [9]:
import sqlite3
import pandas as pd
import matplotlib.pyplot as plt

db_path = os.path.join(rundir, "volumes/volumes.db")
with sqlite3.connect(db_path) as conn:
    # All volumes
    df_all = pd.read_sql_query("SELECT area FROM volumes WHERE area IS NOT NULL", conn)
    # Cluster representatives
    df_rep = pd.read_sql_query("SELECT area FROM volumes WHERE is_representative=1 AND area IS NOT NULL", conn)

areas_all = df_all['area']
areas_rep = df_rep['area']

plt.figure(figsize=(8, 5))
plt.hist(areas_all, bins=40, density=True, color='skyblue', alpha=0.5, label='All volumes')
plt.hist(areas_rep, bins=40, density=True, color='salmon', alpha=0.7, label='Cluster representatives')
plt.xlabel('Area')
plt.ylabel('Density')
plt.yscale('log')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.axvline(np.mean(areas_all), color='blue', linestyle='dashed', linewidth=1, label='Mean All Volumes')
plt.axvline(np.mean(areas_rep), color='red', linestyle='dashed', linewidth=1, label='Mean Representatives')
plt.title('Area Distribution: All Volumes vs Cluster Representatives')
plt.legend()
plt.tight_layout()
plt.show()
No description has been provided for this image
In [6]:
with sqlite3.connect(db_path) as conn:
    # All volumes
    df_all = pd.read_sql_query(
        "SELECT volume, no2d FROM volumes WHERE volume IS NOT NULL AND no2d IS NOT NULL", conn
    )
    # Cluster representatives
    df_rep = pd.read_sql_query(
        "SELECT volume, no2d FROM volumes WHERE is_representative=1 AND volume IS NOT NULL AND no2d IS NOT NULL", conn
    )

plt.figure(figsize=(8, 6))
plt.scatter(df_all['volume'], df_all['no2d'], s=2, color='skyblue', alpha=0.3, label='All volumes')
plt.scatter(df_rep['volume'], df_rep['no2d'], s=40, color='red', alpha=0.8, label='Cluster representatives')
plt.xlabel('Volume')
plt.ylabel('no2d')
plt.title('Volume vs no2d: All Volumes and Cluster Representatives')
plt.legend()
Out[6]:
<matplotlib.legend.Legend at 0x7a3c9ab1c810>
/home/ebr/.cache/pypoetry/virtualenvs/release-volume-sampler-LsGdq2JB-py3.11/lib/python3.11/site-packages/IPython/core/events.py:82: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  func(*args, **kwargs)
/home/ebr/.cache/pypoetry/virtualenvs/release-volume-sampler-LsGdq2JB-py3.11/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image
In [ ]:
# Save write notebook to HTML
!jupyter nbconvert --to html /home/ebr/projects/release-volume-sampler/notebooks/preparational.ipynb --output preparational.html --output-dir "$rundir"
In [ ]: