src/HYDROTools/shapesGroups.py

   1
   2 import salome
   3 salome.salome_init()
   4
   5 import  SMESH, SALOMEDS
   6 from salome.smesh import smeshBuilder
   7
   8 import numpy as np
   9 import MEDLoader as ml
  10 import medcoupling as mc
  11
  12 import shapefile
  13 import math
  14 import os
  15
  16
  17 def freeBordersGroup(meshFileIn, meshFileOut=""):
  18     """
  19     In a mesh, create a group of Edges for the borders: domain, internal holes (isles)
  20     parameters:
  21     meshFileIn: full path of the input mesh file, format MED
  22     meshFileOut: full path of the output mesh file, format MED (default="" : when "", output file is suffixed with "_brd.med"
  23     return full path of the output mesh file
  24     """
  25     smesh = smeshBuilder.New()
  26     smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
  27     ([MESH], status) = smesh.CreateMeshesFromMED(meshFileIn)
  28
  29     nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
  30
  31     aCriteria = []
  32     aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
  33     aCriteria.append(aCriterion)
  34     aFilter = smesh.GetFilterFromCriteria(aCriteria)
  35     aFilter.SetMesh(MESH.GetMesh())
  36     FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
  37
  38     a = os.path.splitext(meshFileIn)
  39     smesh.SetName(MESH, os.path.basename(a[0]))
  40
  41     if meshFileOut == "":
  42         a = os.path.splitext(meshFileIn)
  43         smesh.SetName(MESH, os.path.basename(a[0]))
  44         meshFileOut = a[0] + '_brd' + a[1]
  45
  46     try:
  47         MESH.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
  48         pass
  49     except:
  50         print('ExportMED() failed. Invalid file name?')
  51     return meshFileOut
  52
  53
  54 def explodeGroup(grp, grpName):
  55     """
  56     from a group of edges loaded with MEDCoupling, create ordered lists of nodes, one list for each set of connected edges.
  57     parameters:
  58     grp: MEDCoupling object for the group of edges
  59     grpName: name of the group
  60     return:
  61     List of descriptors [(ordered list of nodeIds, name of the list, closed status)]
  62     """
  63     print(" === explodeGroup", grpName)
  64     nbCells=grp.getNumberOfCells()
  65
  66     dicReverse = {} # id node --> id edges
  67     for i in range(nbCells):
  68         nodcell = grp.getNodeIdsOfCell(i)
  69         for j in range(len(nodcell)):
  70             if nodcell[j] in dicReverse:
  71                 dicReverse[nodcell[j]].append(i)
  72             else:
  73                 dicReverse[nodcell[j]] = [i]
  74
  75     nodeChains = []
  76     usedCells = [False] * nbCells
  77     while False in usedCells:
  78         icell = usedCells.index(False)
  79         usedCells[icell] = True
  80         nodcell = grp.getNodeIdsOfCell(icell)
  81         closed = False
  82         chain = [nodcell[0], nodcell[1]]
  83         nextnode = nodcell[1]
  84         prevnode = nodcell[0]
  85         while nextnode in dicReverse:
  86             nextcells = dicReverse[nextnode]
  87             if len(nextcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  88                 closed = False
  89                 nextnode = -1                   # stop
  90             else:
  91                 newcell =False
  92                 for i in range(len(nextcells)):
  93                     ncell = nextcells[i]
  94                     if not usedCells[ncell]:       # the chain of nodes grows
  95                         usedCells[ncell] = True
  96                         newcell = True
  97                         nodcell = grp.getNodeIdsOfCell(ncell)
  98                         if nodcell[0] == nextnode:
  99                             nextnode = nodcell[1]  # forward edge
 100                         else:
 101                             nextnode = nodcell[0]  # reversed edge ?
 102                         chain.append(nextnode)
 103                 if not newcell:                    # end of chain, closed
 104                     closed =True
 105                     nextnode = -1
 106         while prevnode in dicReverse:
 107             prevcells = dicReverse[prevnode]
 108             if len(prevcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
 109                 closed = False
 110                 prevnode = -1                   # stop
 111             else:
 112                 newcell =False
 113                 for i in range(len(prevcells)):
 114                     ncell = prevcells[i]
 115                     if not usedCells[ncell]:       # the chain of nodes grows
 116                         usedCells[ncell] = True
 117                         newcell = True
 118                         nodcell = grp.getNodeIdsOfCell(ncell)
 119                         if nodcell[1] == prevnode:
 120                             prevnode = nodcell[0]  # forward edge
 121                         else:
 122                             prevnode = nodcell[1]  # reversed edge ?
 123                         chain.insert(0, prevnode)
 124                 if not newcell:                    # end of chain, closed
 125                     closed =True
 126                     prevnode = -1
 127
 128         chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
 129         nodeChains.append(chainDesc)
 130         print(chainDesc[1:])
 131     return nodeChains
 132
 133
 134 def writeShapeLines(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
 135     """
 136     from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type line, with one record per list of connected nodes.
 137     parameters:
 138     mcMesh: mesh loaded in memory with MEDLoader
 139     grpName: name associated to the group of edges
 140     nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
 141     outputDirectory: directory for writing the shapefile
 142     offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
 143     offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
 144     """
 145     print(" === writeShapeLines", grpName)
 146     coords = mcMesh.getCoords()
 147     shapeFileName = os.path.join(outputDirectory, grpName)
 148     w = shapefile.Writer(shapeFileName)
 149     w.shapeType = 3
 150     w.field('name', 'C')
 151     for (chain, chainName, closed) in nodeChains:
 152         print("   --- ", chainName)
 153         chaincoords = []
 154         for node in chain:
 155             coord = coords[node].getValues()
 156             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 157             chaincoords.append(coordLb93)
 158         w.line([chaincoords])
 159         w.record(chainName)
 160     w.close()
 161
 162
 163 def writeShapePoints(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
 164     """
 165     from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type multi points, with one record per list of connected nodes.
 166     parameters:
 167     mcMesh: mesh loaded in memory with MEDLoader
 168     grpName: name associated to the group of edges
 169     nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
 170     outputDirectory: directory for writing the shapefile
 171     offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
 172     offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
 173     """
 174     print(" === writeShapePoints", grpName)
 175     coords = mcMesh.getCoords()
 176     shapeFileName = os.path.join(outputDirectory, grpName)
 177     w = shapefile.Writer(shapeFileName + '_pts')
 178     #w.shapeType = 8
 179     w.field('name', 'C')
 180     for (chain, chainName, closed) in nodeChains:
 181         print("   --- ", chainName)
 182         chaincoords = []
 183         for node in chain:
 184             coord = coords[node].getValues()
 185             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 186             chaincoords.append(coordLb93)
 187         w.multipoint(chaincoords)
 188         w.record(chainName)
 189     w.close()
 190
 191
 192 def exploreEdgeGroups(meshFile, outputDirectory="", offsetX=0., offsetY=0.):
 193     """
 194     Find all the groups of edges in a mesh and, for each group, create one shapefile of lines and one of points. The shapefiles are created in the same system of coordinates as the mesh (For instance Lambert 93), but without origin offset (to be ready for a direct load in Qgis)
 195     parameters:
 196     meshFile: full path of the input mesh file, format MED
 197     outputDirectory: directory in which the shapefiles are written (default "", if "" use the directory containing the mesh
 198     offsetX: local X origin of the mesh
 199     offsetY: local Y origin of the mesh
 200     """
 201     print(" === exploreEdgeGroups", meshFile)
 202     if outputDirectory == "":
 203         outputDirectory = os.path.dirname(meshFile)
 204
 205     a = os.path.splitext(meshFile)
 206     prefix = os.path.basename(a[0]) # prefix = file name without extension
 207
 208     mcMesh = ml.MEDFileMesh.New(meshFile)
 209     dim = mcMesh.getSpaceDimension()
 210     d1=-1        # when dimension 2, edges are dim -1
 211     if dim == 3: # when dimension 3, edges are dim -2
 212         d1=-2
 213
 214     grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
 215
 216     groups = [mcMesh.getGroup(d1, name) for name in grp_names] # list of groups in their name order
 217
 218     for (grp, grpName) in zip(groups, grp_names):
 219         fullGrpName = prefix + '_' + grpName
 220         nodeChains = explodeGroup(grp, fullGrpName)
 221         writeShapeLines(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
 222         writeShapePoints(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
 223
 224
 225 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust, splitFreeBorder=True):
 226     """
 227     shapeFileToAdjust must be a closed line or polygon crossing freeBorderShapefile in 2 points.
 228     Find in shapeFileToAdjust and in freeBorderShapefile the two closest corresponding points and move the points in shapeFileToAdjust to correspond to the points found in freeBorderShapefile. Split shapeFileToAdjust in two parts (inside or outside freeBorder). If requested, split freeBorderShapefile in two parts.
 229     parameters:
 230     freeBorderShapefile: a set of free border lines, as generated by the functions freeBordersGroup and exploreEdgeGroups.
 231     shapefileToAdjust: a closed line or polygon, supposed to be drawn in qgis to pass as close as possible to the points to be connected, on the free border.
 232     splitFreeBorder: boolean default True
 233     """
 234     print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
 235
 236     # --- find domain freeBorder:  bounding box englobing all others
 237     #     TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
 238
 239     fb = shapefile.Reader(freeBorderShapefile)
 240     fbShapes = fb.shapes()
 241     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
 242     outerBboxIndex = -1
 243     for i,s in enumerate(fbShapes):
 244         bb = s.bbox
 245         if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
 246             maxbbox = bb
 247             outerBboxIndex = i
 248     fbs = fbShapes[outerBboxIndex] # the domain free border shape
 249
 250     # --- find the intersections of the shapefile to adjust and the domain free border:
 251     #     the closests points (two pairs of points)
 252
 253     sf = shapefile.Reader(shapefileToAdjust)
 254     shapes = sf.shapes()
 255     sfta = sf.shape(0)
 256     pdist =[]
 257     for i,p in enumerate(sfta.points):
 258         dmin = 1.e30
 259         jmin = -1
 260         for j,pfb in enumerate(fbs.points):
 261             d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
 262             if d < dmin:
 263                 dmin = d
 264                 jmin = j
 265         pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
 266     pdist.sort() # the 2 closest points must be set on the mesh freeBorder
 267     print(pdist)
 268     i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
 269     i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
 270     if i1 == pdist[0][2]:
 271         ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
 272         ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
 273     else:
 274         ifb1 = pdist[1][1]
 275         ifb2 = pdist[0][1]
 276
 277     # --- write the adusted shapefile: free border closest points replaced with corresponding points
 278     #     on the free border. two polylines, one inside the domain, one outside
 279
 280     a = os.path.splitext(shapefileToAdjust)
 281     shapefileAdjusted = a[0] + '_adj' + a[1]
 282     chainName = os.path.basename(a[0])
 283
 284     w = shapefile.Writer(shapefileAdjusted)
 285     w.shapeType = 3
 286     w.field('name', 'C')
 287
 288     chaincoords = []
 289     chaincoords.append(fbs.points[ifb1])
 290     for i in range(i1+1, i2):
 291         chaincoords.append(sfta.points[i])
 292     chaincoords.append(fbs.points[ifb2])
 293     w.line([chaincoords])
 294     w.record(chainName + '_0')
 295
 296     chaincoords = []
 297     chaincoords.append(fbs.points[ifb2])
 298     if i2+1 < len(sfta.points):
 299         for i in range(i2+1, len(sfta.points)):
 300             chaincoords.append(sfta.points[i])
 301     for i in range(i1):
 302         chaincoords.append(sfta.points[i])
 303     chaincoords.append(fbs.points[ifb1])
 304     w.line([chaincoords])
 305     w.record(chainName + '_1')
 306     w.close()
 307
 308     if splitFreeBorder:
 309         # write the free border split in two polylines (cut by the adjusted shapefile)
 310
 311         a = os.path.splitext(freeBorderShapefile)
 312         freeBorderSplit = a[0] + '_split' + a[1]
 313         chainName = os.path.basename(a[0])
 314
 315         w = shapefile.Writer(freeBorderSplit)
 316         w.shapeType = 3
 317         w.field('name', 'C')
 318
 319         if (ifb1 > ifb2):
 320             i = ifb1; ifb1 = ifb2; ifb2 = i
 321
 322         chaincoords = []
 323         for i in range(ifb1, ifb2+1):
 324             chaincoords.append(fbs.points[i])
 325         w.line([chaincoords])
 326         w.record(chainName + '_0')
 327
 328         chaincoords = []
 329         for i in range(ifb2, len(fbs.points)):
 330             chaincoords.append(fbs.points[i])
 331         for i in range(ifb1+1):
 332             chaincoords.append(fbs.points[i])
 333         w.line([chaincoords])
 334         w.record(chainName + '_1')
 335         w.close()
 336