diff --git a/openquake/mbt/tools/tr/catalogue_hmtk.py b/openquake/mbt/tools/tr/catalogue_hmtk.py index 4d22741a0..963d04f8d 100644 --- a/openquake/mbt/tools/tr/catalogue_hmtk.py +++ b/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 diff --git a/openquake/sub/get_profiles_from_slab2pt0.py b/openquake/sub/get_profiles_from_slab2pt0.py index 4595511d5..0e2e36036 100644 --- a/openquake/sub/get_profiles_from_slab2pt0.py +++ b/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 diff --git a/openquake/sub/tests/create_profiles_from_slab2pt0_test.py b/openquake/sub/tests/create_profiles_from_slab2pt0_test.py index 3f3460017..fa3fff544 100644 --- a/openquake/sub/tests/create_profiles_from_slab2pt0_test.py +++ b/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 + + diff --git a/openquake/sub/tests/data/cs_000.csv b/openquake/sub/tests/data/cs_000.csv new file mode 100644 index 000000000..bd4bd3946 --- /dev/null +++ b/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 diff --git a/openquake/sub/tests/data/mar_cs_000.csv b/openquake/sub/tests/data/mar_cs_000.csv new file mode 100644 index 000000000..bd4bd3946 --- /dev/null +++ b/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 diff --git a/openquake/sub/tests/data/vert_dist_to_slab_slab.pkl b/openquake/sub/tests/data/vert_dist_to_slab_slab.pkl new file mode 100644 index 000000000..77bde828c Binary files /dev/null and b/openquake/sub/tests/data/vert_dist_to_slab_slab.pkl differ diff --git a/openquake/sub/tests/profile_workflow_classification_test.py b/openquake/sub/tests/profile_workflow_classification_test.py index 7c10f3f5f..779bee97f 100644 --- a/openquake/sub/tests/profile_workflow_classification_test.py +++ b/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")) diff --git a/openquake/sub/utils.py b/openquake/sub/utils.py index 8545db2e1..51ed27a3e 100644 --- a/openquake/sub/utils.py +++ b/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)