src/HYDROTools/changeBathy.py

   1 # -*- coding: utf-8 -*-
   2
   3 # -----------------------------------------------------------------------------
   4
   5 import numpy as np
   6 import MEDLoader as ml
   7 from scipy import spatial
   8
   9
  10
  11 def changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX=0., offsetY=0., zUndef=-100., interpMethod=0, m3d=False, stats=False):
  12     """
  13     ChangeBathy takes a 2D (x,y) mesh with an already existing field of double ('BOTTOM') for altimetry.
  14     The mesh may already have a z coordinate set with the z altimetry.
  15     The new altimetry from the bathymetry file is set on the nodes of the selected group of face.
  16
  17     In:
  18         meshFileIn: med file name produced by SMESH, corresponding to the HYDRO case
  19         bathyFile: bathymetry file
  20         groupName: name of the group of faces to modify.
  21         offsetX, offsetY: coordinates of origin in MESH
  22         zUndef: Z value to use for nodes outside the bathymetry interpolation (should not happen: if so, algoritm problem).
  23                         default value is -100.
  24         interpMethod: integer value
  25                         0 = nearest point on bathymetry
  26                         1 = linear interpolation
  27         m3d: True/False to produce a 3D mesh. Default is False.
  28         stats: compute statistics for SALOME test
  29     Out:
  30         statz: statistique for z
  31                         Key: face group name
  32                         Value: (minz, maxz, meanz, stdz, v05z, v95z)
  33     Out:
  34         meshFileOut : med file with Z value in a field "BOTTOM"
  35                         Option: Z value also in Z coordinate if m3D is true
  36     """
  37
  38     # --- load the bathymetry in a KD tree, working with the same origin coordinates as the mesh
  39
  40     fi=open(bathyFile, 'r')
  41
  42     minx =  1.e30
  43     maxx = -1.e30
  44     miny =  1.e30
  45     maxy = -1.e30
  46     minz =  1.e30
  47     maxz = -1.e30
  48
  49     coords = []
  50     alti = []
  51     for ligne in fi:
  52         vals = ligne.split()
  53         if len(vals) < 3:
  54             continue
  55         x = float(vals[0]) - offsetX
  56         y = float(vals[1]) - offsetY
  57         z = float(vals[2])
  58         coords.append((x, y))
  59         alti.append(z)
  60         minx = min(x, minx)
  61         miny = min(y, miny)
  62         minz = min(z, minz)
  63         maxx = max(x, maxx)
  64         maxy = max(y, maxy)
  65         maxz = max(z, maxz)
  66     fi.close()
  67
  68     print("bornes x", minx, maxx)
  69     print("bornes y", miny, maxy)
  70     print("bornes z", minz, maxz)
  71     print("nombre points:", len(alti))
  72
  73     bathyPts = np.array(coords)
  74     print(bathyPts.shape)
  75
  76     tree = spatial.KDTree(bathyPts)
  77
  78     # --- load the nodes coordinates of the mesh group
  79
  80     mcMesh = ml.MEDFileMesh.New(meshFileIn)
  81     dim = mcMesh.getMeshDimension()
  82     d0 = 0        # when dimension 2, faces are dim
  83     if dim == 3:  # when dimension 3, edges are dim -1
  84         d0=-1
  85
  86     groupNames = mcMesh.getGroupsOnSpecifiedLev(d0) #names of face groups
  87
  88     if groupName not in groupNames:
  89         print("there is no group of faces named %s in the mesh" %groupName)
  90         #return
  91
  92     group = mcMesh.getGroup(d0, groupName)
  93
  94     # nodes ids in use: number of nodes and order of the mesh, nodes not in the group marked -1,
  95     # nodes indices = 0 to nbNodeInGroup, not their indice in mesh!
  96
  97     AllNodeIds = group.getNodeIdsInUse()[0].toNumPyArray()
  98     print("AllNodeIds.size:", AllNodeIds.size)
  99     cond = AllNodeIds != -1
 100     #nodeIds = np.extract(cond, AllNodeIds)
 101     nids = [i for i in range(AllNodeIds.size) if int(AllNodeIds[i]) >=0] # indices in mesh
 102     nodeIds = np.array(nids)
 103     print("nodeIds.size:", nodeIds.size)
 104     print("nodeIds:", nodeIds)
 105
 106     allCoords = group.getCoords().toNumPyArray()
 107     print("allCoords.shape:",allCoords.shape)
 108     dimCoord = allCoords.shape[1]
 109     if dimCoord == 3:
 110         condition = np.vstack((cond, cond, cond)).transpose()
 111     else:
 112         condition = np.vstack((cond, cond)).transpose()
 113     coords = np.extract(condition, allCoords)
 114     nbRaws = coords.size // dimCoord
 115     coords = coords.reshape(nbRaws, dimCoord)
 116     print("coords.shape:",coords.shape)
 117
 118     if coords.shape[1] == 3:
 119         coordsXY = coords[:, 0:2] # 2 first columns
 120         coordsZ = coords[:, 2] # 3rd column
 121     else:
 122         coordsXY = coords
 123         coordsZ = None
 124
 125     if interpMethod != 0:
 126         print("TODO: the linear interpolation on bathymetry is not yet implemented, using nearest point")
 127
 128     # --- get nearest node Ids for nodes of group
 129
 130     (distances, ids) = tree.query(coordsXY)
 131     print("ids.size", ids.size)
 132     print("ids", ids)
 133
 134     # --- check if there is a 'BOTTOM' Field
 135
 136     try:
 137         bottomField = ml.ReadField(meshFileIn, 'BOTTOM')
 138         bottomArray = bottomField.getArray()
 139     except ml.InterpKernelException:
 140         print("no BOTTOM array")
 141         bottomField = None
 142
 143     if bottomField is None:
 144         bottomField = ml.MEDCouplingFieldDouble(ml.ON_NODES)
 145         bottomField.setName('BOTTOM')
 146         bottomField.setMesh(mcMesh.getMeshAtLevel(0))
 147         npField = np.zeros((mcMesh.getNumberOfNodes()))
 148         bottomArray = ml.DataArrayDouble(npField)
 149
 150     # --- set altitude on field
 151
 152     vz = []
 153     for i,id in enumerate(ids):
 154         #print(i, "---", nodeIds[i], "---", id, alti[id])
 155         j = int(nodeIds[i])
 156         bottomArray[j] = alti[id]
 157         if stats:
 158             vz.append(alti[id])
 159
 160     bottomField.setArray(bottomArray)
 161     bottomField.setTime(0.0, 0, -1)
 162
 163     if (coordsZ is not None) or m3d:
 164         print("TODO: coordsZ")
 165         mlCoords = group.getCoords()
 166         coords3D = ml.DataArrayDouble.Meld([mlCoords[:,0:2], bottomArray])
 167         coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
 168         mcMesh.setCoords(coords3D)
 169
 170     mcMesh.write(meshFileOut, 2)
 171     ml.WriteFieldUsingAlreadyWrittenMesh(meshFileOut, bottomField)
 172
 173     if stats:
 174         statz = {}
 175         minz = np.amin(vz)
 176         maxz = np.amax(vz)
 177         meanz = np.mean(vz)
 178         stdz = np.std(vz)
 179         v05z = np.percentile(vz, 0o5)
 180         v95z = np.percentile(vz, 95)
 181         ligne = ".. Minimum: %f" % minz
 182         ligne += ", maximum: %f" % maxz
 183         ligne += ", mean: %f\n" % meanz
 184         ligne += ".. stdeviation: %f" % stdz
 185         ligne += ", v05z: %f" % v05z
 186         ligne += ", v95z: %f" % v95z
 187         print (ligne)
 188         statz[groupName] = (minz, maxz, meanz, stdz, v05z, v95z)
 189         return statz
 190
 191
 192 if __name__=='__main__':
 193     meshFileIn = 'garonne_1F.med'
 194     meshFileOut = 'garonne_1FF.med'
 195     bathyFile = 'newBathy.xyz'
 196     #groupName = 'litMineur'
 197     groupName = 'riveGauche'
 198     #groupName = 'riveDroite'
 199     offsetX = 430000
 200     offsetY = 6350000
 201     stats = changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX, offsetY, stats=True)
 202
 203