# Export of the calculation case
garonne_1_entry = garonne_1.Export()
+# --- add a new bathymetry for the test changeBathy
+
+import tempfile
+tmpdir = tempfile.mkdtemp()
+print("tmpdir=",tmpdir)
+
+newBathy = os.path.join(tmpdir, 'newBathy.xyz')
+fi=open(os.path.join(HYDRO_SAMPLES, "garonne_point_L93.xyz" ), 'r')
+fo=open(newBathy, 'w')
+for ligne in fi:
+ vals = ligne.split()
+ if len(vals) < 3:
+ continue
+ x = float(vals[0])
+ y = float(vals[1])
+ z = float(vals[2]) + 50
+ l = "%12.3f %12.3f %12.3f\n" % (x, y, z)
+ fo.write(l)
+fi.close()
+fo.close()
+
#----------------------
# --- Geometry
#----------------------
edges_riveGaucheIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveGauche[i]) for i in range(len(edges_riveGauche)) ]
edges_riveDroiteIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveDroite[i]) for i in range(len(edges_riveDroite)) ]
-print("edges_litMineurIds", edges_litMineurIds)
+print("edges_litMineurIds", edges_litMineurIds)
print("edges_riveGaucheIds", edges_riveGaucheIds)
print("edges_riveDroiteIds", edges_riveDroiteIds)
'litMineur': (2.06, 25.41, 13.93, 4.33, 8.47, 21.78)}
controlStatZ(statz, refstatz)
-# --- add a field on nodes of type double with z values, named "BOTTOM"
-#createZfield2(fichierMaillage)
+# --- interpolation on a new bathymetry for 'riveGauche', without using HYDRO (no Calculation Case)
+
+meshFileIn = os.path.splitext(fichierMaillage)[0] + 'F.med'
+meshFileOut = os.path.splitext(meshFileIn)[0] + 'F.med'
+# --- Z interpolation on the bathymety/altimetry on the mesh nodes
+
+from salome.hydrotools.changeBathy import changeBathy
+statz = changeBathy(meshFileIn, meshFileOut, newBathy, 'riveGauche', 430000, 6350000, stats=True)
+print(statz)
+refstatz = {'riveGauche': (57.72, 77.14, 71.33, 3.36, 62.92, 75.0)}
+controlStatZ(statz, refstatz)
--- /dev/null
+# -*- coding: utf-8 -*-
+
+# -----------------------------------------------------------------------------
+
+import numpy as np
+import MEDLoader as ml
+from scipy import spatial
+
+
+
+def changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX=0., offsetY=0., zUndef=-100., interpMethod=0, m3d=False, stats=False):
+ """
+ ChangeBathy takes a 2D (x,y) mesh with an already existing field of double ('BOTTOM') for altimetry.
+ The mesh may already have a z coordinate set with the z altimetry.
+ The new altimetry from the bathymetry file is set on the nodes of the selected group of face.
+
+ In:
+ meshFileIn: med file name produced by SMESH, corresponding to the HYDRO case
+ bathyFile: bathymetry file
+ groupName: name of the group of faces to modify.
+ offsetX, offsetY: coordinates of origin in MESH
+ zUndef: Z value to use for nodes outside the bathymetry interpolation (should not happen: if so, algoritm problem).
+ default value is -100.
+ interpMethod: integer value
+ 0 = nearest point on bathymetry
+ 1 = linear interpolation
+ m3d: True/False to produce a 3D mesh. Default is False.
+ stats: compute statistics for SALOME test
+ Out:
+ statz: statistique for z
+ Key: face group name
+ Value: (minz, maxz, meanz, stdz, v05z, v95z)
+ Out:
+ meshFileOut : med file with Z value in a field "BOTTOM"
+ Option: Z value also in Z coordinate if m3D is true
+ """
+
+ # --- load the bathymetry in a KD tree, working with the same origin coordinates as the mesh
+
+ fi=open(bathyFile, 'r')
+
+ minx = 1.e30
+ maxx = -1.e30
+ miny = 1.e30
+ maxy = -1.e30
+ minz = 1.e30
+ maxz = -1.e30
+
+ coords = []
+ alti = []
+ for ligne in fi:
+ vals = ligne.split()
+ if len(vals) < 3:
+ continue
+ x = float(vals[0]) - offsetX
+ y = float(vals[1]) - offsetY
+ z = float(vals[2])
+ coords.append((x, y))
+ alti.append(z)
+ minx = min(x, minx)
+ miny = min(y, miny)
+ minz = min(z, minz)
+ maxx = max(x, maxx)
+ maxy = max(y, maxy)
+ maxz = max(z, maxz)
+ fi.close()
+
+ print("bornes x", minx, maxx)
+ print("bornes y", miny, maxy)
+ print("bornes z", minz, maxz)
+ print("nombre points:", len(alti))
+
+ bathyPts = np.array(coords)
+ print(bathyPts.shape)
+
+ tree = spatial.KDTree(bathyPts)
+
+ # --- load the nodes coordinates of the mesh group
+
+ mcMesh = ml.MEDFileMesh.New(meshFileIn)
+ dim = mcMesh.getMeshDimension()
+ d0 = 0 # when dimension 2, faces are dim
+ if dim == 3: # when dimension 3, edges are dim -1
+ d0=-1
+
+ groupNames = mcMesh.getGroupsOnSpecifiedLev(d0) #names of face groups
+
+ if groupName not in groupNames:
+ print("there is no group of faces named %s in the mesh" %groupName)
+ #return
+
+ group = mcMesh.getGroup(d0, groupName)
+
+ # nodes ids in use: number of nodes and order of the mesh, nodes not in the group marked -1,
+ # nodes indices = 0 to nbNodeInGroup, not their indice in mesh!
+
+ AllNodeIds = group.getNodeIdsInUse()[0].toNumPyArray()
+ print("AllNodeIds.size:", AllNodeIds.size)
+ cond = AllNodeIds != -1
+ #nodeIds = np.extract(cond, AllNodeIds)
+ nids = [i for i in range(AllNodeIds.size) if int(AllNodeIds[i]) >=0] # indices in mesh
+ nodeIds = np.array(nids)
+ print("nodeIds.size:", nodeIds.size)
+ print("nodeIds:", nodeIds)
+
+ allCoords = group.getCoords().toNumPyArray()
+ print("allCoords.shape:",allCoords.shape)
+ dimCoord = allCoords.shape[1]
+ if dimCoord == 3:
+ condition = np.vstack((cond, cond, cond)).transpose()
+ else:
+ condition = np.vstack((cond, cond)).transpose()
+ coords = np.extract(condition, allCoords)
+ nbRaws = coords.size // dimCoord
+ coords = coords.reshape(nbRaws, dimCoord)
+ print("coords.shape:",coords.shape)
+
+ if coords.shape[1] == 3:
+ coordsXY = coords[:, 0:2] # 2 first columns
+ coordsZ = coords[:, 2] # 3rd column
+ else:
+ coordsXY = coords
+ coordsZ = None
+
+ if interpMethod != 0:
+ print("TODO: the linear interpolation on bathymetry is not yet implemented, using nearest point")
+
+ # --- get nearest node Ids for nodes of group
+
+ (distances, ids) = tree.query(coordsXY)
+ print("ids.size", ids.size)
+ print("ids", ids)
+
+ # --- check if there is a 'BOTTOM' Field
+
+ try:
+ bottomField = ml.ReadField(meshFileIn, 'BOTTOM')
+ bottomArray = bottomField.getArray()
+ except ml.InterpKernelException:
+ print("no BOTTOM array")
+ bottomField = None
+
+ if bottomField is None:
+ bottomField = ml.MEDCouplingFieldDouble(ml.ON_NODES)
+ bottomField.setName('BOTTOM')
+ bottomField.setMesh(mcMesh.getMeshAtLevel(0))
+ npField = np.zeros((mcMesh.getNumberOfNodes()))
+ bottomArray = ml.DataArrayDouble(npField)
+
+ # --- set altitude on field
+
+ vz = []
+ for i,id in enumerate(ids):
+ #print(i, "---", nodeIds[i], "---", id, alti[id])
+ j = int(nodeIds[i])
+ bottomArray[j] = alti[id]
+ if stats:
+ vz.append(alti[id])
+
+ bottomField.setArray(bottomArray)
+ bottomField.setTime(0.0, 0, -1)
+
+ if (coordsZ is not None) or m3d:
+ print("TODO: coordsZ")
+ mlCoords = group.getCoords()
+ coords3D = ml.DataArrayDouble.Meld([mlCoords[:,0:2], bottomArray])
+ coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
+ mcMesh.setCoords(coords3D)
+
+ mcMesh.write(meshFileOut, 2)
+ ml.WriteFieldUsingAlreadyWrittenMesh(meshFileOut, bottomField)
+
+ if stats:
+ statz = {}
+ minz = np.amin(vz)
+ maxz = np.amax(vz)
+ meanz = np.mean(vz)
+ stdz = np.std(vz)
+ v05z = np.percentile(vz, 0o5)
+ v95z = np.percentile(vz, 95)
+ ligne = ".. Minimum: %f" % minz
+ ligne += ", maximum: %f" % maxz
+ ligne += ", mean: %f\n" % meanz
+ ligne += ".. stdeviation: %f" % stdz
+ ligne += ", v05z: %f" % v05z
+ ligne += ", v95z: %f" % v95z
+ print (ligne)
+ statz[groupName] = (minz, maxz, meanz, stdz, v05z, v95z)
+ return statz
+
+
+if __name__=='__main__':
+ meshFileIn = 'garonne_1F.med'
+ meshFileOut = 'garonne_1FF.med'
+ bathyFile = 'newBathy.xyz'
+ #groupName = 'litMineur'
+ groupName = 'riveGauche'
+ #groupName = 'riveDroite'
+ offsetX = 430000
+ offsetY = 6350000
+ stats = changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX, offsetY, stats=True)
+
+
