#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Meshing tools
"""
import os
# from scipy.spatial.distance import cdist
# from scipy.spatial import KDTree
import numpy as np
import pyvista as pv
pv.set_plot_theme("document")
#pv.set_jupyter_backend('static')
import pandas as pd
from pyCATHY.plotters import cathy_plots as cplt
from scipy.spatial import KDTree
try:
import pygimli as pg
import pygimli.meshtools as mt
except ImportError:
pygimli = None
import matplotlib.pylab as plt
def CATHY_2_Simpeg(
mesh_CATHY, ERT_meta_dict, scalar="saturation", show=False, **kwargs
):
pass
def CATHY_2_pg(mesh_CATHY, ERT_meta_dict, scalar="saturation", show=False, **kwargs):
"""
Interpolate CATHY mesh attribute to pygimli mesh.
Add a new [`scalar`] attribute to the pygimli mesh (create a new mesh)
.. Note:
Need to flip axis because convention for CATHY and pygimli are different
Parameters
----------
mesh_CATHY : pvmesh
CATHY mesh to transform to pygimli.
ERT_meta_dict :dict
Dictionnary containing ERT metadata (mesh, format, ..).
scalar : str, optional
scalar attribute to interpolate. The default is 'saturation'.
show : bool, optional
show the result of the interpolation using pyvista. The default is False.
**kwargs : TYPE
path : path of the mesh to overwrite
Returns
-------
mesh_new_attr : TYPE
DESCRIPTION.
scalar_new : TYPE
DESCRIPTION.
"""
if type(ERT_meta_dict["forward_mesh_vtk_file"]) is str:
mesh_OUT = pv.read(ERT_meta_dict["forward_mesh_vtk_file"])
# flip y and z axis as CATHY and pg have different convention for axis
# ------------------------------------------------------------------------
in_nodes_mod = np.array(mesh_CATHY.points)
# in_nodes_mod_pg = np.array(mesh_OUT.points)
# mesh_nodes_modif = None
if 'mesh_nodes_modif' in ERT_meta_dict:
print('mesh transformation before interpolation')
in_nodes_mod_m = ERT_meta_dict['mesh_nodes_modif']
else:
print('no mesh transformation before interpolation')
in_nodes_mod_m = in_nodes_mod[:, :]
path = os.getcwd()
if "path" in kwargs:
path = kwargs["path"]
meshCATHY_tmp = mesh_CATHY.copy()
data_OUT, warm_0 = trace_mesh(
meshCATHY_tmp,
mesh_OUT,
scalar=scalar,
threshold=1e-1,
in_nodes_mod=in_nodes_mod_m
)
if len(warm_0) > 0:
print(warm_0)
print('add new attribute to pg mesh')
scalar_new = scalar + "_nearIntrp2_pg_msh"
if "time" in kwargs:
time = kwargs["time"]
mesh_new_attr, name_new_attr = add_attribute_2mesh(
data_OUT, mesh_OUT, scalar_new, overwrite=True, time=time, path=path
)
else:
mesh_new_attr, name_new_attr = add_attribute_2mesh(
data_OUT, mesh_OUT, scalar_new, overwrite=True, path=path
)
print('end of CATHY_2_pg')
show = False
if show:
p = pv.Plotter(window_size=[1024 * 3, 768 * 2], off_screen=True,) #notebook=True)
p.add_mesh(mesh_new_attr, scalars=scalar_new)
_ = p.add_bounding_box(line_width=5, color="black")
cpos = p.show(True)
p.save_graphic(
'test21.svg', title="",
# raster=True,
# painter=True
)
p = pv.Plotter(window_size=[1024 * 3, 768 * 2], off_screen=True,) #notebook=True)
p.add_mesh(mesh_CATHY, scalars=scalar)
_ = p.add_bounding_box(line_width=5, color="black")
cpos = p.show(True)
return mesh_new_attr, scalar_new
def CATHY_2_Resipy(mesh_CATHY, mesh_Resipy, scalar="saturation", show=False, **kwargs):
# flip y and z axis as CATHY and Resipy have different convention for axis
# ------------------------------------------------------------------------
in_nodes_mod = np.array(mesh_CATHY.points)
in_nodes_mod[:, 2] = -np.flipud(in_nodes_mod[:, 2])
in_nodes_mod[:, 1] = -np.flipud(in_nodes_mod[:, 1])
# check with a plot positon of the nodes for both meshes
# ------------------------------------------------------------------------
# p = pv.Plotter(window_size=[1024*3, 768*2], notebook=True)
# p.add_mesh(mesh_CATHY)
# _ = p.add_points(np.array(mesh_CATHY.points), render_points_as_spheres=True,
# color='red', point_size=20)
# _ = p.add_points(in_nodes_mod, render_points_as_spheres=True,
# color='blue', point_size=20)
# _ = p.show_bounds(grid='front', all_edges=True, font_size=50)
# cpos = p.show(True)
path = os.getcwd()
if "path" in kwargs:
path = kwargs["path"]
data_OUT = trace_mesh(
mesh_CATHY,
mesh_Resipy,
scalar=scalar,
threshold=1e-1,
# threshold=1e2,
in_nodes_mod=in_nodes_mod,
)
# print(len(data_OUT))
scalar_new = scalar + "_nearIntrp2Resipymsh"
print('Add new attribute to mesh')
# get_array(mesh, name, preference='cell'
if "time" in kwargs:
time = kwargs["time"]
mesh_new_attr, name_new_attr = add_attribute_2mesh(
data_OUT, mesh_Resipy, scalar_new, overwrite=True, time=time, path=path
)
else:
mesh_new_attr, name_new_attr = add_attribute_2mesh(
data_OUT, mesh_Resipy, scalar_new, overwrite=True, path=path
)
if show == True:
p = pv.Plotter(window_size=[1024 * 3, 768 * 2], off_screen=True,) #notebook=True)
p.add_mesh(mesh_new_attr, scalars=scalar_new)
_ = p.add_bounding_box(line_width=5, color="black")
cpos = p.show(True)
p.save_graphic(
'test21.svg', title="",
raster=True,
painter=True
)
p = pv.Plotter(window_size=[1024 * 3, 768 * 2], off_screen=True,) #notebook=True)
p.add_mesh(mesh_CATHY, scalars=scalar)
_ = p.add_bounding_box(line_width=5, color="black")
cpos = p.show(True)
# if type(meshERT) is str:
# meshERTpv = pv.read(meshERT)
# if savefig == True:
# plotter = pv.Plotter(notebook=True)
# _ = plotter.add_mesh(mesh_new_attr,show_edges=True)
# plotter.view_xz(negative=False)
# plotter.show_grid()
p.save_graphic(
'test.svg', title="",
raster=True,
painter=True
)
return mesh_new_attr, scalar_new
[docs]
def trace_mesh(meshIN, meshOUT, scalar, threshold=1e-1, **kwargs):
"""
Trace meshIN on meshOUT using nearest neigbour interpolation
Parameters
----------
meshIN : TYPE
DESCRIPTION.
meshOUT : TYPE
DESCRIPTION.
threshold : TYPE, optional
DESCRIPTION. The default is 1e-1.
Returns
-------
out_data : TYPE
DESCRIPTION.
"""
in_nodes_mod = kwargs["in_nodes_mod"]
meshIN.set_active_scalars(scalar)
# meshIN.points = in_nodes_mod
meshOUT.points = in_nodes_mod
# len(in_nodes_mod)
# print(meshIN.points)
# set_interpolation_radius()
rd = max(np.diff(meshIN.points[:, 0])) / 1
meshOUT_interp = meshOUT.interpolate(meshIN,
radius=rd,
pass_point_data=True
)
# plot_2d_interpolation_quality(meshIN,scalar,meshOUT,meshOUT_interp)
result = meshOUT_interp.point_data_to_cell_data()
out_data = result[scalar]
# out_data = np.where(out_data == 0, 1e-3, out_data)
out_data = np.where(out_data == 0, np.min(out_data)+1e-3, out_data)
warm_0 = ""
if len(np.where(out_data == 1e-3)) > 0:
warm_0 = f"interpolation created 0 values - replacing them by min value {np.min(out_data)+1e-3} of input CATHY predicted ER mesh"
warm_0 = f"min {np.min(out_data)}, max {np.min(out_data)}, median {np.median(out_data)} "
return out_data, warm_0
[docs]
def set_interpolation_radius():
# rd= min([abs(min(np.diff(meshIN.points[:,0]))),
# abs(min(np.diff(meshIN.points[:,1]))),
# abs(min(np.diff(meshIN.points[:,2])))
# ]
# )
pass
[docs]
def plot_2d_interpolation_quality(meshIN, scalar, meshOUT, result):
# fig = plt.figure()
# ax1 = plt.subplot(131)
# # print(max(meshIN[scalar]))
# # print(min(meshIN[scalar]))
# # meshIN.points[:,0].min()
# # meshOUT.points[:,0].min()
# # meshOUT.points[:,0].max()
# # meshOUT.points[:,1].min()
# # meshOUT.points[:,1].max()
# cm = plt.cm.get_cmap('RdYlBu')
# # result = meshOUT.interpolate(meshIN, radius=rd, pass_point_data=True)
# sc = ax1.scatter(meshIN.points[:,0],meshIN.points[:,1],c=meshIN[scalar],label='meshIN[scalar]',
# s=35, cmap=cm)
# plt.colorbar(sc)
cm = plt.cm.get_cmap("RdYlBu")
# # fig = plt.figure()
fig, axs = plt.subplots(2)
# ax2 = plt.subplot(121)
sc = axs[0].scatter(
meshIN.points[:, 0],
meshIN.points[:, 1],
c=meshIN[scalar],
label="meshIN[scalar]",
s=35,
cmap=cm,
)
# ,vmin=min(meshIN[scalar]), vmax=max(meshIN[scalar])
axs[0].set_xlim([min(meshIN.points[:, 0]), max(meshIN.points[:, 0])])
axs[0].set_ylim([min(meshIN.points[:, 1]), max(meshIN.points[:, 1])])
# fig = plt.figure()
# ax3 = plt.subplot(122)
sc = axs[1].scatter(
meshOUT.points[:, 0],
meshOUT.points[:, 1],
c=result[scalar],
label="result[scalar]",
s=35,
cmap=cm,
)
axs[1].set_xlim([min(meshIN.points[:, 0]), max(meshIN.points[:, 0])])
axs[1].set_ylim([min(meshIN.points[:, 1]), max(meshIN.points[:, 1])])
import matplotlib.ticker as ticker
def fmt(x, pos):
a, b = "{:.2e}".format(x).split("e")
b = int(b)
return r"${} \times 10^{{{}}}$".format(a, b)
plt.colorbar(sc, format=ticker.FuncFormatter(fmt))
# plt.tight_layout()
def uniquify(path):
filename, extension = os.path.splitext(path)
counter = 1
while os.path.exists(path):
path = filename + str(counter) + extension
counter += 1
return path
savedir = os.getcwd()
savename_test = os.path.join(savedir, "interpolation_q.png")
savename = uniquify(savename_test)
# print(savename)
plt.savefig(savename)
# plt.show()
pass
# def find_nearest_cellcenter(node_coord,meshIN_nodes_coords,threshold=1e-1,
# **kwargs):
# '''
# Find nearest mesh node between two meshes
# Parameters
# ----------
# node_coord : np.array
# meshIN_nodes_coords : np.array
# threshold : float
# if distance > threshold --> closest = nan
# Returns
# -------
# closest_idx : list
# Node indice in the mesh_node.
# closest : list
# Node coordinate in the mesh_node.
# '''
# closest_idx = []
# closest = []
# # for i, nc in enumerate(cell_coords):
# # euclidean distance
# d = ( (meshIN_nodes_coords[:,0] - node_coord[0]) ** 2 +
# (meshIN_nodes_coords[:,1] - node_coord[1]) ** 2 +
# (meshIN_nodes_coords[:,2] - node_coord[2]) ** 2
# ) ** 0.5
# closest_idx.append(np.argmin(d))
# closest.append(np.vstack(meshIN_nodes_coords[closest_idx,:]))
# if d[np.argmin(d)]>threshold:
# closest = 'nan'
# return closest_idx, closest
[docs]
def find_nearest_node(node_coord, meshIN_nodes_coords, threshold=1e-1, **kwargs):
"""
Find nearest mesh node between two meshes
Parameters
----------
node_coord : np.array
meshIN_nodes_coords : np.array
threshold : float
if distance > threshold --> closest = nan
Returns
-------
closest_idx : list
Node indice in the mesh_node.
closest : list
Node coordinate in the mesh_node.
"""
closest_idx = []
closest = []
# for i, nc in enumerate(cell_coords):
# euclidean distance
d = (
(meshIN_nodes_coords[:, 0] - node_coord[0]) ** 2
+ (meshIN_nodes_coords[:, 1] - node_coord[1]) ** 2
+ (meshIN_nodes_coords[:, 2] - node_coord[2]) ** 2
) ** 0.5
closest_idx.append(np.argmin(d))
closest.append(np.vstack(meshIN_nodes_coords[closest_idx, :]))
if d[np.argmin(d)] > threshold:
closest = "nan"
return closest_idx, closest
[docs]
def add_attribute_2mesh(
data, mesh, name="ER_pred", overwrite=True, saveMesh=False, **kwargs
):
"""
add a new mesh attribute to a vtk file
Parameters
----------
data : TYPE
DESCRIPTION.
mesh : TYPE
DESCRIPTION.
name : TYPE
DESCRIPTION.
overwrite : TYPE, optional
DESCRIPTION. The default is True.
**kwargs : TYPE
DESCRIPTION.
Returns
-------
None.
"""
# for k in kwargs:
# print(k)
# print(mesh)
if type(mesh) is str:
mesh = pv.read(mesh)
try:
mesh.point_data[name] = data
except:
mesh.cell_data[name] = data
meshname = name + ".vtk"
saveMesh = True
if saveMesh:
path = os.getcwd()
if "path" in kwargs:
path = kwargs["path"]
if "time" in kwargs:
time = kwargs["time"]
meshname = name + str(time) + ".vtk"
mesh.save(path + meshname, binary=False)
else:
mesh.save(path + meshname, binary=False)
print(path + meshname)
# if overwrite==True:
# mesh.save(full_path)
return mesh, name
[docs]
def trace_mesh_pg(meshIN, meshOUT, method="spline", **kwargs):
"""
Interpolate CATHY mesh (structured) into pygimli mesh (structured) using pygimli meshtools
# Specify interpolation method 'linear, 'spline', 'harmonic'
"""
meshIN = pg.load(meshIN)
meshOUT = pg.load(meshOUT)
out_data = pg.interpolate(meshIN["ER_converted_CATHYmsh*"], meshOUT, method=method)
return out_data
#%%
def map_layers_2_DEM(layers, DEM, zone, dem_parameters):
ltop, lbot = get_layer_depths(dem_parameters)
zone3d_layers_top, zone3d_layer_bot = get_zone3d_layer_depths(zone, dem_parameters)
dem_mat3d_layers_top = [DEM - zz for zz in zone3d_layers_top]
dem_mat3d_layers_bot = [DEM - zz for zz in zone3d_layer_bot]
zone3d_topflag = []
for li in range(dem_parameters["nstr"]):
zone3d_topflag_li = np.ones(np.shape(dem_mat3d_layers_top[0]))
bool_top_lli = []
for ll in layers.keys():
layers_adj_top = DEM - abs(layers[ll][0])
layers_adj_bot = DEM - abs(layers[ll][1])
# differences between top of the layer i of the mesh and top of the desired layer
# -------------------------------------------------------------------------------
diff_top = dem_mat3d_layers_top[li] - layers_adj_top
cond1 = diff_top <= 1e-2
# -------------------------------------------------------------------------------
diff_bot = dem_mat3d_layers_bot[li] - layers_adj_bot
cond2 = diff_bot <= 1e-2
# print('lmeshi'+str(li),
# ll,layers[ll],np.mean(diff_top),np.mean(diff_bot),cond1[0][0],cond2[0][0])
# if li<dem_parameters["nstr"]-1:
# cond2 = abs(diff)<= abs(dem_mat3d_layers_top[li+1] - dem_mat3d_layers_top[li])
# else:
# cond2 = np.ones(np.shape(diff),dtype=bool)
bool_top_lli = cond1 & cond2
zone3d_topflag_li[bool_top_lli] = ll
zone3d_topflag.append(zone3d_topflag_li)
zone3d_topflag = np.array(zone3d_topflag)
return zone3d_topflag
def get_zone3d_layer_depths(zone_raster, dem_parameters):
'''
Return a 3d matrice with dimension [Number of layers, X cells, Y cells]
Parameters
----------
zone_raster : TYPE
DESCRIPTION.
dem_parameters : TYPE
DESCRIPTION.
Returns
-------
TYPE
DESCRIPTION.
'''
zone3d_top = []
zone3d_bot = []
zone3d_raster = np.array([zone_raster] * dem_parameters["nstr"])
for li in range(dem_parameters["nstr"]):
top, bot = get_layer_depth(dem_parameters, li)
zone3d_top.append(zone3d_raster[li] * top)
zone3d_bot.append(zone3d_raster[li] * bot)
return np.array(zone3d_top), np.array(zone3d_bot)
def get_layer_depths(dem_parameters):
layers_top = []
layers_bottom = []
for li in range(dem_parameters["nstr"]):
layeri_top, layeri_bottom = get_layer_depth(dem_parameters, li)
layers_top.append(layeri_top)
layers_bottom.append(layeri_bottom)
return layers_top, layers_bottom
def get_layer_depth(dem_parameters, li):
if type(dem_parameters["zratio(i),i=1,nstr"]) != list:
dempar = dem_parameters["zratio(i),i=1,nstr"].split("\t")
dempar_ratio = [float(d) for d in dempar]
else:
dempar_ratio = [float(d) for d in dem_parameters["zratio(i),i=1,nstr"]]
if li == 0:
layeri_top = 0
else:
layeri_top = np.cumsum(dempar_ratio[0:li])[-1] * (dem_parameters["base"])
if (li + 1) < len(dempar_ratio):
layeri_bottom = np.cumsum(dempar_ratio[0 : li + 1])[-1] * (
dem_parameters["base"]
)
else:
layeri_bottom = dem_parameters["base"]
return layeri_top, layeri_bottom
def zone3d(zone, dem_parameters):
# define zone in 3dimension - duplicate number of layer (nstr) times
# ---------------------------------------------------------------
zones3d = [zone] * dem_parameters["nstr"]
# fig, axs = plt.subplots(
# int(dem_parameters["nstr"] / 2) + 1,
# 2,
# sharex=True,
# sharey=(True),
# constrained_layout=False,
# )
# plt.tight_layout()
# axs = axs.ravel()
# for li in range(dem_parameters["nstr"]):
# layeri_top, layeri_bottom = get_layer_depth(dem_parameters, li)
# layer_str = (
# "[" + str(f"{layeri_top:.2E}") + "-" + str(f"{layeri_bottom:.2E}") + "]"
# )
# pmesh = cplt.show_raster(
# zones3d[li], prop=layer_str, cmap="jet", ax=axs[li], vmin=0, vmax=1
# )
# plt.colorbar(pmesh, ax=axs[:], location="right", shrink=0.6, cmap="jet")
# plt.tight_layout()
return zones3d
def create_layers_inzones3d(simu, zones3d, layers_names, layers_depths=[[0, 1e99]]):
# Loop over layers and zones and change flag if depth conditions is not respected
# ----------------------------------------------------------------------------
zones3d_layered = np.ones(np.shape(zones3d))
for li in range(simu.dem_parameters["nstr"]):
layeri_top, layeri_bottom = get_layer_depth(simu, li)
# layer_str = '[' + str(layeri_top) + '-' + str(layeri_bottom) + ']'
for zi in range(len(layers_names)):
print("layers %i analysis" % zi)
# print(layers_depths[zi][0])
# print(layeri_top)
# print(layers_depths[zi][1])
# print(layeri_bottom)
zi = 0
# print(layers_depths[zi][0]<=layeri_top)
# if depth of zone i is sup to mesh layers depth
# ---------------------------------------------------------------------
# if (layers_depths[zi][0]>=layeri_top) & (layers_depths[zi][1]<layeri_bottom):
if (layers_depths[zi][0] <= layeri_top) & (
layers_depths[zi][1] >= layeri_bottom
):
# if (depths_ordered[zi][1]>=layeri_bottom):
print("conds ok --> zi:" + str(zi))
print(layers_depths[zi])
print(layeri_top, layeri_bottom)
zones3d_layered[li][zones3d[li] == zi + 1] = zi + 2
print(
"replacing "
+ str(np.sum(zones3d[li] == zi + 1))
+ " values by"
+ str(zi + 1)
)
if 10.5 * layers_depths[zi][1] < layeri_bottom:
raise ValueError(
"Required layer is finer than mesh layers - refine mesh"
)
else:
print("conds not ok")
print(layers_depths[zi])
print(layeri_top, layeri_bottom)
return zones3d_layered
def plot_zones3d_layered(simu, zones3d_layered):
fig, axs = plt.subplots(
int(simu.dem_parameters["nstr"] / 2) + 1,
2,
sharex=False,
sharey=(False),
constrained_layout=False,
)
axs = axs.flat
for li in range(simu.dem_parameters["nstr"]):
layeri_top, layeri_bottom = get_layer_depth(simu, li)
layer_str = "[" + str(layeri_top) + "-" + str(layeri_bottom) + "]"
pmesh = cplt.show_raster(
zones3d_layered[li],
prop=layer_str, # , cmap='jet',
ax=axs[li],
vmin=0,
vmax=2,
)
plt.colorbar(pmesh, ax=axs[:], location="right", shrink=0.6, cmap="jet")
plt.tight_layout()
def het_soil_layers_mapping_generic(
simu, propertie_names, SPP, layers_names, layers_depths
):
# extend to 3d the zone raster file
# -----------------------------------
zones3d = zone3d(simu)
# insert layers flag into the 3d the zone raster file
# -----------------------------------
zones3d_layered = create_layers_inzones3d(
simu, zones3d, layers_names, layers_depths
)
# plot 3d zones files layered
# ------------------------------------------
plot_zones3d_layered(simu, zones3d_layered)
# np.shape(zones3d_layered)
# np.shape(zones3d_axis_swap)
# Loop over soil properties names
# -------------------------------------------------
index_raster = np.arange(0, simu.hapin["N"] * simu.hapin["M"])
zones3d_axis_swap = np.swapaxes(zones3d_layered, 0, 2)
prop_df = [] # properties dataframe
layers_id = ["L" + str(i + 1) for i in range(simu.dem_parameters["nstr"])]
layers_id = np.flipud(layers_id)
for i, p in enumerate(propertie_names):
p_df = np.zeros(np.shape(zones3d_axis_swap))
# Loop over soil layers and assign value of soil properties
# -------------------------------------------------
for k, lname in enumerate(layers_names):
p_df[zones3d_axis_swap == k + 1] = SPP[k][i]
df_tmp = pd.DataFrame(
np.vstack(p_df),
columns=layers_id,
index=index_raster,
)
prop_df.append(df_tmp)
SoilPhysProp_df_het_layers_p = pd.concat(prop_df, axis=1, keys=propertie_names)
SoilPhysProp_df_het_layers_p.index.name = "id raster"
SoilPhysProp_df_het_layers_p.columns.names = ["soilp", "layerid"]
SPP_map_dict = {}
for p in propertie_names:
SPP_map_dict[p] = []
for li in range(simu.dem_parameters["nstr"]):
v = SoilPhysProp_df_het_layers_p.iloc[0][p].loc["L" + str(li + 1)]
SPP_map_dict[p].append(v)
return SoilPhysProp_df_het_layers_p, SPP_map_dict
def _subplot_cellsMarkerpts(mesh_pv_attributes, xyz_layers0, xyzlayers1):
pl = pv.Plotter(shape=(1, 2))
# pl.add_mesh(mesh_pv_attributes,
# show_edges=True,
# )
pl.show_grid()
actor = pl.add_points(
xyz_layers0[:, :-1],
point_size=10,
scalars=xyz_layers0[:, -1],
)
pl.subplot(0, 1)
pl.show_grid()
actor = pl.add_points(
xyzlayers1[:, :-1],
point_size=10,
scalars=xyzlayers1[:, -1],
)
pl.show()
def _plot_cellsMarkerpts(mesh_pv_attributes, xyz_layers):
pl = pv.Plotter()
pl.add_mesh(mesh_pv_attributes, show_edges=True, opacity=0.4)
pl.show_grid()
actor = pl.add_points(
xyz_layers[:, :-1],
point_size=10,
scalars=xyz_layers[:, -1],
)
pl.set_scale(zscale=5)
pl.show()
def _build_xyz_marker_mat_squareDEM(
zones3d_layered,
mesh_pv_attributes,
hapin,
dem_parameters,
workdir,
project_name,
to_nodes,
):
''' DEPRECATED FOR _find_nearest_point2DEM'''
# top, bot = get_layer_depth(dem_parameters,li)
# get layers properties
# ----------------------------------------------------------------
dempar = dem_parameters["zratio(i),i=1,nstr"].split("\t")
dempar_ratio = [float(d) for d in dempar]
layeri_top = []
layeri_top = [
np.cumsum(dempar_ratio[0 : li + 1])[-1] * (dem_parameters["base"])
for li in range(dem_parameters["nstr"])
]
layeri_top.insert(0, 0)
layeri_center = [
(layeri_top[lti + 1] + layeri_top[lti]) / 2
for lti in range(len(layeri_top) - 1)
]
# get dem coordinates
# ----------------------------------------------------------------
y, x = cplt.get_dem_coords(dem, hapin=hapin)
# build dem coordinates of the size of the DEM
# ----------------------------------------------------------------
xn = [np.ones(len(y)) * xu for xu in np.unique(x)]
xn = np.hstack(xn)
yn = np.tile(y, len(x))
xyz = np.c_[
xn - dem_parameters["delta_x"] / 2,
yn - dem_parameters["delta_y"] / 2,
np.hstack(dem.T),
]
xyz_celldown = np.c_[
xn - dem_parameters["delta_x"] / 4,
yn - dem_parameters["delta_y"] / 4,
np.hstack(dem.T),
]
xyz_cellup = np.c_[
xn + dem_parameters["delta_x"] / 4,
yn + dem_parameters["delta_y"] / 4,
np.hstack(dem.T),
]
xyz_cells = np.r_[xyz_cellup, xyz_celldown]
# len(xyz_cellup)
xyz_nodes = np.c_[xn, yn, np.hstack(dem.T)]
# meshpts = mesh_pv_attributes.points
# cellpts = mesh_pv_attributes.cell_centers().points
xyz_layers_nodes = []
for i, li in enumerate(layeri_top[:-1]):
marker_zone = np.ravel(zones3d_layered[i])
xyz_layers_nodes.append(
np.c_[xyz_nodes[:, 0:2], xyz_nodes[:, 2] - li, marker_zone]
)
xyz_layers = []
for i, li in enumerate(layeri_center[:]):
marker_zone = np.ravel(zones3d_layered[i])
xyz_layers.append(np.c_[xyz[:, 0:2], xyz[:, 2] - li, marker_zone])
xyz_layers_cells_top = []
for i, li in enumerate(layeri_top[:-1]):
marker_zone = np.hstack([np.ravel(zones3d_layered[i])] * 2)
xyz_layers_cells_top.append(
np.c_[xyz_cells[:, 0:2], xyz_cells[:, 2] - li, marker_zone]
)
xyz_layers_cells_cent = []
for i, li in enumerate(layeri_center[:]):
marker_zone = np.hstack([np.ravel(zones3d_layered[i])] * 2)
xyz_layers_cells_cent.append(
np.c_[xyz_cells[:, 0:2], xyz_cells[:, 2] - li, marker_zone]
)
xyz_layers = np.vstack(xyz_layers)
xyz_layers_nodes = np.vstack(xyz_layers_nodes)
xyz_layers_cells = np.vstack(np.r_[xyz_layers_cells_top,
xyz_layers_cells_cent]
)
_plot_cellsMarkerpts(mesh_pv_attributes, xyz_layers, workdir, project_name)
_find_nearest_point2DEM(
to_nodes,
mesh_pv_attributes,
xyz_layers_cells,
xyz_layers_nodes,
)
def _find_nearest_point2DEM(
to_nodes,
mesh_pv_attributes,
xyz_layers_cells=[],
xyz_layers_nodes=[],
saveMeshPath=None,
):
if to_nodes:
# loop over mesh cell centers and find nearest point to dem
# ----------------------------------------------------------------
node_markers = []
for nmesh in mesh_pv_attributes.points:
# euclidean distance
d = (
(xyz_layers_nodes[:, 0] - nmesh[0]) ** 2
+ (xyz_layers_nodes[:, 1] - nmesh[1]) ** 2
+ (xyz_layers_nodes[:, 2] - nmesh[2]) ** 2
)** 0.5
node_markers.append(xyz_layers_nodes[np.argmin(d), 3])
# add data to the mesh
# ----------------------------------------------------------------
mesh_pv_attributes["node_markers_old"] = node_markers
else:
# Get the mesh cell centers
cell_centers = mesh_pv_attributes.cell_centers().points
# Create a KDTree from xyz_layers_cells
tree = KDTree(xyz_layers_cells[:, :3])
# Find the nearest neighbors for each cell center
distances, indices = tree.query(cell_centers, k=1)
# Get the cell_markers values based on the nearest neighbors
cell_markers = xyz_layers_cells[indices, 3]
# Add the data to the mesh
mesh_pv_attributes["cell_markers"] = cell_markers
nodepv = mesh_pv_attributes.cell_data_to_point_data()
noderounded = [np.ceil(npv) for npv in nodepv.point_data['cell_markers']]
mesh_pv_attributes['node_markers_new'] = noderounded
if saveMeshPath is not None:
mesh_pv_attributes.save(saveMeshPath,binary=False)
def add_markers2mesh(
zones3d_layered,
dem,
mesh_pv_attributes,
dem_parameters,
hapin,
grid3d,
to_nodes=False,
show=False,
saveMeshPath=None,
):
# Reduce all to 1D
# ------------------------------------------------------------------
zones3d_layers_array = np.ravel(zones3d_layered)
x, y = cplt.get_dem_coords(dem,hapin=hapin)
dem_flip = dem
# Get layer top and bottom
# ------------------------------------------------------------------
zone3d_top = []
zone3d_bot = []
for li in range(dem_parameters["nstr"]):
top, bot = get_layer_depth(dem_parameters, li)
zone3d_top.append(np.ones(np.shape(zones3d_layered[li]))*top)
zone3d_bot.append(np.ones(np.shape(zones3d_layered[li]))*bot)
zone3d_top = np.array(zone3d_top)
zone3d_bot = np.array(zone3d_bot)
if to_nodes:
dem_mat3d_layers = [dem_flip - zz for zz in zone3d_top-(zone3d_top-zone3d_bot)/2]
else:
dem_mat3d_layers = [dem_flip - zz for zz in zone3d_top-(zone3d_top-zone3d_bot)/2]
# Create a regular mesh from the DEM x and y coordinates and elevation
# ------------------------------------------------------------------
xgrid, ygrid = np.meshgrid(x, y)
grid_coords_dem = np.array(
[
np.ravel(xgrid),
np.ravel(ygrid),
]
).T
# Reduce all to 1D
# ------------------------------------------------------------------
dem_mat_stk = np.ravel(dem_mat3d_layers)
# dem_mat_stk = np.fliplr(dem_mat_stk)
grid_coords_stk_rep = np.vstack(np.array([grid_coords_dem] * dem_parameters["nstr"]))
zones3d_col_stk = np.ravel(zones3d_layered)
xyz_layers = np.c_[grid_coords_stk_rep, dem_mat_stk, zones3d_col_stk]
# Plot to check position of points VS mesh
# ------------------------------------------------------------------
if show:
_plot_cellsMarkerpts(mesh_pv_attributes,
xyz_layers,
)
#%%
# print('_find_nearest_point2DEM')
# Assign marker to mesh and overwrite it
# ------------------------------------------------------------------
_find_nearest_point2DEM(
to_nodes,
mesh_pv_attributes,
xyz_layers,
xyz_layers,
saveMeshPath
)