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