WIP: Classifier updates for 2026

openquake/mbt/tools/tr/catalogue_hmtk.py

@@ -61,9 +61,10 @@ def get_distances_from_surface(catalogue, surface):
             mesh = Mesh(catalogue.data['longitude'][i:upp],
                         catalogue.data['latitude'][i:upp],
                         catalogue.data['depth'][i:upp])
-            tmp = surface.get_min_distance(mesh)
+            tmp = np.float64(surface.get_min_distance(mesh), 5)
             dsts[i:upp] = tmp
             i = upp 
+    dsts = dsts.round(5)
     return dsts
 
 

openquake/sub/get_profiles_from_slab2pt0.py

@@ -113,7 +113,7 @@ def get_initial_traces(box, boy, dip_dir, spacing):
     # reverse profiles to satisfy right hand rule
     profiles.reverse()
 
-    return np.array(profiles), max_length
+    return np.asarray(profiles), max_length
 
 
 def aa_get_initial_traces(bb, dip_dir, spacing):
@@ -145,24 +145,24 @@ def aa_get_initial_traces(bb, dip_dir, spacing):
 
     points = g.fwd_intermediate(bb[0], bb[2], az12_low, num_samples_low + 1,
                                 spacing_hor, initial_idx=0, terminus_idx=0)
-    profiles = _get_profiles(np.array(points.lons), np.array(points.lats),
+    profiles = _get_profiles(np.asarray(points.lons), np.asarray(points.lats),
                              dip_dir, max_length)
 
     points = g.fwd_intermediate(bb[1], bb[3], az21_right,
                                 num_samples_right + 1, spacing_ver,
                                 initial_idx=0, terminus_idx=0)
-    tmp_profiles = _get_profiles(np.array(points.lons), np.array(points.lats),
+    tmp_profiles = _get_profiles(np.asarray(points.lons), np.asarray(points.lats),
                                  dip_dir, max_length)
     profiles.extend(tmp_profiles)
 
-    return np.array(profiles), max_length
+    return np.asarray(profiles), max_length
 
 
 def _get_profiles(lons, lats, dip_dir, max_length):
     tmp_profiles = []
     for xco, yco in zip(lons, lats):
         tmp = point_at(xco, yco, dip_dir, max_length)
-        arr = np.array([[xco, yco], [tmp[0], tmp[1]]])
+        arr = np.asarray([[xco, yco], [tmp[0], tmp[1]]])
         tmp_profiles.append(arr)
     return tmp_profiles
 
@@ -244,7 +244,7 @@ def tmp_get_initial_traces(bb, dip_dir, spacing):
 
         # Update the profile list
         profiles.extend(tmp_profiles)
-        profiles = np.array(profiles)
+        profiles = np.asarray(profiles)
         mask = profiles[:, :, 0] > 180
         profiles[mask, 0] = profiles[mask, 0] - 360
         return profiles, distance
@@ -300,12 +300,12 @@ def get_profiles_geojson(geojson: str, fname_dep: str, spacing: float,
 
     # Reading file with depth values
     f_dep = netCDF4.Dataset(fname_dep)
-    depths = np.array(f_dep.variables['z'])
+    depths = np.asarray(f_dep.variables['z'])
     mask = np.where(np.isfinite(depths))
 
     # Mesh
-    x = np.array(f_dep.variables['x'])
-    y = np.array(f_dep.variables['y'])
+    x = np.asarray(f_dep.variables['x'])
+    y = np.asarray(f_dep.variables['y'])
     xx, yy = np.meshgrid(x, y)
 
     css = []
@@ -318,7 +318,7 @@ def get_profiles_geojson(geojson: str, fname_dep: str, spacing: float,
     max_lo = -180.0
     max_la = -90.0
     for index, row in gdf.iterrows():
-        coo = np.array(row.coords)
+        coo = np.asarray(row.coords)
         min_lo = np.min([min_lo, np.min(coo[:, 0])])
         min_la = np.min([min_la, np.min(coo[:, 1])])
         max_lo = np.max([max_lo, np.max(coo[:, 0])])
@@ -348,7 +348,7 @@ def get_profiles_geojson(geojson: str, fname_dep: str, spacing: float,
     # Coords
     tmp = zip(xx.flatten(), yy.flatten(), depths.flatten())
     depths = [[x, y, z] for x, y, z in tmp]
-    depths = np.array(depths)
+    depths = np.asarray(depths)
     mask = depths[:, 0] > 180
     depths[mask, 0] = depths[mask, 0] - 360
     milo = np.min(depths[:, 0])
@@ -362,7 +362,7 @@ def get_profiles_geojson(geojson: str, fname_dep: str, spacing: float,
 
 
     if len(str(fname_fig)) > 0:
-        bb = np.array([min_lo, min_la, max_lo, max_la])
+        bb = np.asarray([min_lo, min_la, max_lo, max_la])
         dlt = 5.0
         reg = [bb[0] - dlt, bb[2] + dlt, bb[1] - dlt, bb[3] + dlt]
         clo = np.mean([bb[0], bb[1]])
@@ -428,7 +428,7 @@ def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
 
     # Reading file with strike values
     f_strike = netCDF4.Dataset(fname_str)
-    strikes = np.array(f_strike.variables['z'])
+    strikes = np.asarray(f_strike.variables['z'])
     mask = np.where(np.isfinite(strikes))
     strikes = strikes[mask]
 
@@ -437,8 +437,8 @@ def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
     dip_dir = (strike_dir + 90) % 360
 
     # Mesh
-    x = np.array(f_strike.variables['x'])
-    y = np.array(f_strike.variables['y'])
+    x = np.asarray(f_strike.variables['x'])
+    y = np.asarray(f_strike.variables['y'])
     xx, yy = np.meshgrid(x, y)
 
     # Compute the initial bounding box
@@ -471,7 +471,7 @@ def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
 
     # Reading file with depth values
     f_dep = netCDF4.Dataset(fname_dep)
-    depths = np.array(f_dep.variables['z'])
+    depths = np.asarray(f_dep.variables['z'])
     mask = np.where(np.isfinite(depths))
 
     # Filter
@@ -482,7 +482,7 @@ def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
     # Coords
     tmp = zip(xx.flatten(), yy.flatten(), depths.flatten())
     depths = [[x, y, z] for x, y, z in tmp]
-    depths = np.array(depths)
+    depths = np.asarray(depths)
     mask = depths[:, 0] > 180
     depths[mask, 0] = depths[mask, 0] - 360
 

openquake/sub/tests/create_profiles_from_slab2pt0_test.py

@@ -27,6 +27,9 @@
 import pathlib
 import unittest
 import numpy as np
+import pandas as pd
+import tempfile
+import os
 from openquake.sub.get_profiles_from_slab2pt0 import (get_profiles,
                                                       get_bounding_box,
                                                       get_initial_traces,
@@ -200,4 +203,16 @@ def test_geojson_mariana(self):
         fname = HERE / 'data' / 'izu_slab2_css.geojson'
         fname_slab = HERE / 'data' / 'izu_slab2_dep_02.24.18.grd'
 
-        get_profiles_geojson(fname, fname_slab, spacing=10.)
+        test = get_profiles_geojson(fname, fname_slab, spacing=10.)
+
+        top_level = tempfile.mkdtemp()
+        profiles = tempfile.mkdtemp(dir=top_level)
+
+        test.write_profiles(profiles)
+
+        df1 = pd.read_csv(os.path.join(profiles,'cs_000.csv'))
+        df2 = pd.read_csv(HERE / "data" / "mar_cs_000.csv")
+        res = df1.compare(df2)
+        assert len(res) == 0
+
+

openquake/sub/tests/data/cs_000.csv

@@ -0,0 +1,49 @@
+1.435446094052026638e+02 2.401308591980515672e+01 1.132732022828012930e+01
+1.434954976799403141e+02 2.400434697468638490e+01 1.190561589960553590e+01
+1.434463926316922198e+02 2.399559238759614033e+01 1.254143236962554298e+01
+1.433972942687359762e+02 2.398682216709252302e+01 1.323666306944010174e+01
+1.433482025993247078e+02 2.397803632174514021e+01 1.399429457634280638e+01
+1.432991176316870394e+02 2.396923486013503179e+01 1.482556792962662051e+01
+1.432500393740271534e+02 2.396041779085468448e+01 1.574420432198734154e+01
+1.432009678345243913e+02 2.395158512250799276e+01 1.676243620142352242e+01
+1.431519030213338510e+02 2.394273686371022336e+01 1.789757172483989933e+01
+1.431028449425858184e+02 2.393387302308800102e+01 1.917136028048950180e+01
+1.430537936063860514e+02 2.392499360927928720e+01 2.058329299147684566e+01
+1.430047490208158649e+02 2.391609863093334454e+01 2.216642013631120278e+01
+1.429557111939317338e+02 2.390718809671070133e+01 2.394920448041169792e+01
+1.429066801337656614e+02 2.389826201528315508e+01 2.592329399299599046e+01
+1.428576558483248675e+02 2.388932039533371210e+01 2.814328151808559753e+01
+1.428086383455921293e+02 2.388036324555658751e+01 3.066574464814294032e+01
+1.427596276335253265e+02 2.387139057465717684e+01 3.346200916374149870e+01
+1.427106237200578960e+02 2.386240239135200270e+01 3.655092701965818947e+01
+1.426616266130984343e+02 2.385339870436873255e+01 4.002451856989834766e+01
+1.426126363205308678e+02 2.384437952244611836e+01 4.387326164097233061e+01
+1.425636528502144529e+02 2.383534485433397165e+01 4.814472404329134037e+01
+1.425146762099837758e+02 2.382629470879317424e+01 5.280716198140613926e+01
+1.424657064076485540e+02 2.381722909459560711e+01 5.795058157771575225e+01
+1.424167434509938630e+02 2.380814802052413981e+01 6.356865348581769126e+01
+1.423677873477800517e+02 2.379905149537262687e+01 6.962007400241589039e+01
+1.423188381057426852e+02 2.378993952794584033e+01 7.616785668034182777e+01
+1.422698957325925448e+02 2.378081212705949099e+01 8.316454153878673594e+01
+1.422209602360155714e+02 2.377166930154015745e+01 9.059205530690070418e+01
+1.421720316236730639e+02 2.376251106022529314e+01 9.842075024588027077e+01
+1.421231099032013674e+02 2.375333741196317661e+01 1.066362707307830959e+02
+1.420741950822122419e+02 2.374414836561291509e+01 1.152026983293106071e+02
+1.420252871682923512e+02 2.373494393004436986e+01 1.240593132900936837e+02
+1.419763861690037743e+02 2.372572411413819893e+01 1.331799377969740021e+02
+1.419274920918835505e+02 2.371648892678576459e+01 1.425402296540906946e+02
+1.418786049444440778e+02 2.370723837688915125e+01 1.521112089443612376e+02
+1.418297247341726859e+02 2.369797247336111568e+01 1.618230276431121410e+02
+1.417808514685321200e+02 2.368869122512506920e+01 1.716798226686515250e+02
+1.417319851549599434e+02 2.367939464111505998e+01 1.817318634691140460e+02
+1.416831258008690213e+02 2.367008273027573040e+01 1.918661884548160401e+02
+1.416342734136473780e+02 2.366075550156229923e+01 2.022366833295921822e+02
+1.415854280006580268e+02 2.365141296394055104e+01 2.127530311619178747e+02
+1.415365895692390552e+02 2.364205512638677931e+01 2.234935955372048113e+02
+1.414877581267038806e+02 2.363268199788778290e+01 2.345162610869746800e+02
+1.414389336803406820e+02 2.362329358744083763e+01 2.458748628952209856e+02
+1.413901162374128546e+02 2.361388990405365007e+01 2.576196954502950121e+02
+1.413413058051589530e+02 2.360447095674437890e+01 2.698029403007594169e+02
+1.412925023907924356e+02 2.359503675454154603e+01 2.825116878153254447e+02
+1.412437060015019483e+02 2.358558730648406154e+01 2.957778391302533691e+02
+1.411949166444511263e+02 2.357612262162117744e+01 3.096842848306729366e+02

openquake/sub/tests/data/mar_cs_000.csv

@@ -0,0 +1,49 @@
+1.435446094052026638e+02 2.401308591980515672e+01 1.132732022828012930e+01
+1.434954976799403141e+02 2.400434697468638490e+01 1.190561589960553590e+01
+1.434463926316922198e+02 2.399559238759614033e+01 1.254143236962554298e+01
+1.433972942687359762e+02 2.398682216709252302e+01 1.323666306944010174e+01
+1.433482025993247078e+02 2.397803632174514021e+01 1.399429457634280638e+01
+1.432991176316870394e+02 2.396923486013503179e+01 1.482556792962662051e+01
+1.432500393740271534e+02 2.396041779085468448e+01 1.574420432198734154e+01
+1.432009678345243913e+02 2.395158512250799276e+01 1.676243620142352242e+01
+1.431519030213338510e+02 2.394273686371022336e+01 1.789757172483989933e+01
+1.431028449425858184e+02 2.393387302308800102e+01 1.917136028048950180e+01
+1.430537936063860514e+02 2.392499360927928720e+01 2.058329299147684566e+01
+1.430047490208158649e+02 2.391609863093334454e+01 2.216642013631120278e+01
+1.429557111939317338e+02 2.390718809671070133e+01 2.394920448041169792e+01
+1.429066801337656614e+02 2.389826201528315508e+01 2.592329399299599046e+01
+1.428576558483248675e+02 2.388932039533371210e+01 2.814328151808559753e+01
+1.428086383455921293e+02 2.388036324555658751e+01 3.066574464814294032e+01
+1.427596276335253265e+02 2.387139057465717684e+01 3.346200916374149870e+01
+1.427106237200578960e+02 2.386240239135200270e+01 3.655092701965818947e+01
+1.426616266130984343e+02 2.385339870436873255e+01 4.002451856989834766e+01
+1.426126363205308678e+02 2.384437952244611836e+01 4.387326164097233061e+01
+1.425636528502144529e+02 2.383534485433397165e+01 4.814472404329134037e+01
+1.425146762099837758e+02 2.382629470879317424e+01 5.280716198140613926e+01
+1.424657064076485540e+02 2.381722909459560711e+01 5.795058157771575225e+01
+1.424167434509938630e+02 2.380814802052413981e+01 6.356865348581769126e+01
+1.423677873477800517e+02 2.379905149537262687e+01 6.962007400241589039e+01
+1.423188381057426852e+02 2.378993952794584033e+01 7.616785668034182777e+01
+1.422698957325925448e+02 2.378081212705949099e+01 8.316454153878673594e+01
+1.422209602360155714e+02 2.377166930154015745e+01 9.059205530690070418e+01
+1.421720316236730639e+02 2.376251106022529314e+01 9.842075024588027077e+01
+1.421231099032013674e+02 2.375333741196317661e+01 1.066362707307830959e+02
+1.420741950822122419e+02 2.374414836561291509e+01 1.152026983293106071e+02
+1.420252871682923512e+02 2.373494393004436986e+01 1.240593132900936837e+02
+1.419763861690037743e+02 2.372572411413819893e+01 1.331799377969740021e+02
+1.419274920918835505e+02 2.371648892678576459e+01 1.425402296540906946e+02
+1.418786049444440778e+02 2.370723837688915125e+01 1.521112089443612376e+02
+1.418297247341726859e+02 2.369797247336111568e+01 1.618230276431121410e+02
+1.417808514685321200e+02 2.368869122512506920e+01 1.716798226686515250e+02
+1.417319851549599434e+02 2.367939464111505998e+01 1.817318634691140460e+02
+1.416831258008690213e+02 2.367008273027573040e+01 1.918661884548160401e+02
+1.416342734136473780e+02 2.366075550156229923e+01 2.022366833295921822e+02
+1.415854280006580268e+02 2.365141296394055104e+01 2.127530311619178747e+02
+1.415365895692390552e+02 2.364205512638677931e+01 2.234935955372048113e+02
+1.414877581267038806e+02 2.363268199788778290e+01 2.345162610869746800e+02
+1.414389336803406820e+02 2.362329358744083763e+01 2.458748628952209856e+02
+1.413901162374128546e+02 2.361388990405365007e+01 2.576196954502950121e+02
+1.413413058051589530e+02 2.360447095674437890e+01 2.698029403007594169e+02
+1.412925023907924356e+02 2.359503675454154603e+01 2.825116878153254447e+02
+1.412437060015019483e+02 2.358558730648406154e+01 2.957778391302533691e+02
+1.411949166444511263e+02 2.357612262162117744e+01 3.096842848306729366e+02

openquake/sub/tests/profile_workflow_classification_test.py

@@ -13,12 +13,14 @@
 
 
 def test_geojson_classification_workflow():
-    raise unittest.SkipTest
     fname_geojson = HERE / "data" / "izu_slab2_css.geojson"
     fname_slab = HERE / "data" / "izu_slab2_dep_02.24.18.grd"
     classification_ini = HERE / "data" / "mariana_classification.ini"
 
     ## these go in a tmp folder...
+    #top_level =  HERE / "data" / "test_folder"
+    #profiles = HERE / "data" / "test_folder" / "profiles"
+    #sfc_out = HERE / "data" / "test_folder" / "sfc_out"
     top_level = tempfile.mkdtemp()
     profiles = tempfile.mkdtemp(dir=top_level)
     sfc_out = tempfile.mkdtemp(dir=top_level)
@@ -43,13 +45,23 @@ def test_geojson_classification_workflow():
     with open(classification_ini, "w") as configfile:
         config.write(configfile)
 
-    max_sampl_dist = 25
-    upper_depth_int = 0
-    lower_depth_int = 50
-    lower_depth_slab = 300
+    max_sampl_dist = 25.0
+    upper_depth_int = 0.0
+    lower_depth_int = 50.0
+    lower_depth_slab = 300.0
 
     slb = get_profiles_geojson(fname_geojson, fname_slab, spacing=10.0)
     slb.write_profiles(profiles)
+
+    # Compare slab profiles
+    import pandas as pd
+
+    df1 = pd.read_csv(os.path.join(profiles, "cs_000.csv"))
+    df2 = pd.read_csv(HERE / "data" / "cs_000.csv")
+
+    # Compare DataFrames
+    res = df1.compare(df2)
+    assert len(res) == 0
 
     create_2pt5_model(profiles, sfc_out, max_sampl_dist)
 
@@ -68,6 +80,13 @@ def test_geojson_classification_workflow():
     )
 
     classify.classify(classification_ini, True)
+
+    # Test the distances are the same 
+    dists = pd.read_pickle(os.path.join(top_level, "distances", "vert_dist_to_slab_slab.pkl"))
+    test_dists = pd.read_pickle(HERE / "data" / "vert_dist_to_slab_slab.pkl")
+    assert len(dists) == len(test_dists)
+    for a, b in zip(dists, test_dists):
+        assert a == b
 
     # Test the resulting classification is as expected
     classified = pd.read_csv(os.path.join(top_level, "classified_earthquakes.csv"))

openquake/sub/utils.py

@@ -306,6 +306,7 @@ def _check_edges(edges):
     # Fit the plane
     tmp = np.vstack((x.flatten(), y.flatten(), pnts[:, 2].flatten())).T
     _, ppar = plane_fit(tmp)
+    pparray = np.array([ppar[0], ppar[1], 0])  
 
     # Analyse the edges
     chks = []
@@ -319,8 +320,9 @@ def _check_edges(edges):
 
         # Check the edge direction Vs the perpendicular to the plane
         idx = [0, -1]
-        edgv = np.array([np.diff(ex[idx])[0], np.diff(ey[idx])[0]])
-        chks.append(np.sign(np.cross(ppar[:2], edgv)))
+        edgv = np.array([np.diff(ex[idx])[0], np.diff(ey[idx])[0], 0])
+        chks.append(np.sign(np.cross(pparray, edgv)[2]))
+
 
     return np.array(chks)