src/HYDROTools/hydroGeoMeshUtils.py

   1 import sys
   2 import os
   3 import salome
   4
   5 salome.salome_init()
   6
   7 from HYDROPy import *
   8 from PyQt5.QtCore import *
   9 from PyQt5.QtGui import *
  10
  11 from random import randint
  12
  13 def getChildrenInStudy(obj):
  14     """
  15     Given an object published in SALOME study (for instance a GEOM object), retreive its children.
  16     parameters:
  17     obj: a SALOME object published in SALOME study
  18     return:
  19     a dictionary [name] --> object
  20     """
  21     SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(obj))
  22     childIterator = salome.myStudy.NewChildIterator(SO)
  23     children = {}
  24     while childIterator.More():
  25         childItem = childIterator.Value()
  26         print("item", childItem)
  27         itemName = childItem.GetName()
  28         itemID = childItem.GetID()
  29         itemObj = salome.IDToObject(str(itemID))
  30         children[itemName] = itemObj
  31         childIterator.Next()
  32     return children
  33
  34 def loadImage(document, imageFile, imageName, displayLevel):
  35     """
  36     load an image (photo, map...) in order to georeference it and use it as background for drawing shapes
  37     parameters:
  38     document: current HYDROData document
  39     imageFile: full path of the image file to load
  40     imageName: name to give to the object loaded in the HYDRO document
  41     displayLevel: integer >=0, the higher levels are displayed on top
  42     return:
  43     the object image loaded
  44     """
  45     image = document.CreateObject(KIND_IMAGE)
  46     image.SetName(imageName)
  47     image.SetZLevel(displayLevel)
  48     if not(image.LoadImage(imageFile)):
  49         raise ValueError('problem while loading image %s'%imageFile)
  50     image.Update()
  51     return image
  52
  53 def GeolocaliseImageCoords(image, a, b, c ,d):
  54     """
  55     Geolocalize a loaded image with to sets of points: pixels coordinates <--> absolutes coordinates in the geographic reference system.
  56     parameters:
  57     image: loaded image to localize
  58     a: first pixel coordinates QPoint(x, y) (integers)
  59     b: second pixel coordinates QPoint(x, y) (integers)
  60     c: first geographic point QPointF(X, Y) (floats)
  61     d: second geographic point QPointF(X, Y) (floats)
  62     """
  63     image.SetLocalPoints(a, b)
  64     image.SetGlobalPoints(1, c, d)
  65     image.Update()
  66
  67 def GeolocaliseImageReference(image, imageRef, a, b, c ,d):
  68     """
  69     Geolocalize an image using another one already localized, with two pair of corresponding pixels on the images
  70     parameters:
  71     image: loaded image to localize
  72     imageRef: loaded image used as a reference
  73     a, b : pixel coordinates on the image to localize QPoint(x, y) (integers)
  74     c, d : pixel coordinates on the reference image QPoint(x, y) (integers)
  75     """
  76     image.SetLocalPoints(a, b)
  77     image.SetGlobalPoints(3, c, d)
  78     image.SetTrsfReferenceImage(imageRef)
  79     image.Update()
  80
  81 def importPolylines(document, shapeFile, iSpline, displayLevel):
  82     """
  83     Import a QGis shape (*.shp) containing one or several 2D polylines or polygons.
  84     The polyline is either named from its name in the shape database if existing, or using the shapefile basename
  85     with a suffix '_PolyXY'.
  86     In case of several polylines, name are suffixed in sequence with '_n'.
  87     parameters:
  88     document: current HYDROData document
  89     shapeFile: full path of the shapefile (*.shp)
  90     isSpline: boolean to convert the polylines in  splines
  91     displayLevel: integer >=0, the higher levels are displayed on top
  92     return:
  93     list of shape objects
  94     """
  95     entities = HYDROData_PolylineXY.ImportShapesFromFile(shapeFile)
  96     shapeType = 0 # polyline
  97     if iSpline:
  98         shapeType = 1 # polyline
  99     shapes = []
 100     for shape in entities:
 101         print(shape.GetName())
 102         if shape.GetKind() == KIND_POLYLINEXY:
 103             print("polyline %s"%shape.GetName())
 104             for i in range(shape.NbSections()):
 105                 shape.SetSectionType(i, shapeType)
 106                 shape.Update()
 107                 shape.SetZLevel(displayLevel)
 108             shapes.append(shape)
 109     return shapes
 110
 111 def splitShapeTool(document, aShape, toolShape, tolerance=1.e-2):
 112     """
 113     Split the shape (2D polyline or spline) by the shape tool.
 114     return a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 115     parameters;
 116     document: current HYDROData document
 117     aShape: a loaded shape (2D polyline or spline)
 118     toolShape: a loaded shape (2D polyline or spline)
 119     tolerance: used by the algorithm to detect intersections, default value 1.e-2
 120     return:
 121     a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 122     """
 123     name = aShape.GetName()
 124     existingNbSplit = {}
 125     index = 1
 126     found = True
 127     while found:
 128         nameIndex = name + "_%d"%(index)
 129         shape = document.FindObjectByName(nameIndex)
 130         if shape is None:
 131             found = False
 132         else:
 133             index = index + 1
 134     existingNbSplit[name] = index
 135     print(existingNbSplit)
 136
 137     op = HYDROData_PolylineOperator()
 138     op.SplitTool(document, aShape, toolShape, tolerance)
 139
 140     listSplit = []
 141     index = existingNbSplit[name]
 142     found = True
 143     while found:
 144         nameIndex = name + "_%d"%(index)
 145         shape = document.FindObjectByName(nameIndex)
 146         if shape is None:
 147             found = False
 148         else:
 149             print("found %s"%nameIndex)
 150             index = index + 1
 151             listSplit.append(shape)
 152     return listSplit
 153
 154 def splitShapesAll(document, shapeslist, precision=1.E-2):
 155     """
 156     Split all the shapes (2D polylines or splines) in the list by all the other shapes.
 157     return a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 158     parameters;
 159     document: current HYDROData document
 160     shapeslist: a list of loaded shapes (2D polylines or splines)
 161     precision: used by the algorithm to detect intersections, default value 1.e-2
 162     return:
 163     a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 164     """
 165     names = [s.GetName() for s in shapeslist]
 166
 167     existingNbSplit = {}
 168     for name in names:
 169         index = 1
 170         found = True
 171         while found:
 172             nameIndex = name + "_%d"%(index)
 173             shape = document.FindObjectByName(nameIndex)
 174             if shape is None:
 175                 found = False
 176             else:
 177                 index = index + 1
 178         existingNbSplit[name] = index
 179     print(existingNbSplit)
 180
 181     op = HYDROData_PolylineOperator()
 182     op.SplitAll(document, shapeslist, precision)
 183
 184     listSplit = []
 185     for name in names:
 186         index = existingNbSplit[name]
 187         found = True
 188         while found:
 189             nameIndex = name + "_%d"%(index)
 190             shape = document.FindObjectByName(nameIndex)
 191             if shape is None:
 192                 found = False
 193             else:
 194                 print("found %s"%nameIndex)
 195                 index = index + 1
 196                 listSplit.append(shape)
 197     return listSplit
 198
 199
 200 def importBathymetry(document, bathyFile):
 201     """
 202     import a bathymetry file (*.xyz or *.asc)
 203     parameters:
 204     document: current HYDROData document
 205     bathyFile: full path of the bathymetry (*.xyz or *.shp)
 206     return:
 207     the bathymetry object
 208     """
 209     bathy = document.CreateObject(KIND_BATHYMETRY)
 210     a = os.path.splitext(bathyFile)
 211     bathyName = os.path.basename(a[0])
 212     bathy.SetName(bathyName)
 213     bathy.SetAltitudesInverted(0)
 214     if not(bathy.ImportFromFile(bathyFile)):
 215         raise ValueError('problem while loading bathymetry')
 216     bathy.Update()
 217     return bathy
 218
 219 def createImmersibleZone(document, imzName, polyLine, bathy, isImmersible, displayLevel):
 220     """
 221     Create an immersible or unsubmersible zone with a closed polyline and a bathymetry.
 222     parameters:
 223     document: current HYDROData document
 224     imzName: the name to give to the zone object
 225     polyline: the closed 2D polyline or spline delimiting the zone
 226     bathy: the bathymetry to associate to the zone
 227     isImmersible: boolean
 228     displayLevel: integer >=0, the higher levels are displayed on top
 229     return:
 230     the zone created
 231     """
 232     imz = document.CreateObject(KIND_IMMERSIBLE_ZONE)
 233     imz.SetName(imzName)
 234     imz.SetZLevel(displayLevel)
 235     if bathy is not None:
 236         imz.SetAltitudeObject(bathy)
 237     imz.SetPolyline(polyLine)
 238     imz.SetIsSubmersible(isImmersible)
 239     imz.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 240     imz.Update()
 241     return imz
 242
 243 def mergePolylines(document, polyName, polyLines, isConnectedBySegment = False, tolerance = 1.E-3):
 244     """
 245     Regroup several 2D polylines or spline in one, with several sections
 246     parameters:
 247     document: current HYDROData document
 248     polyName: name to give to the new object
 249     polyLines: a list of loaded shapes (2D polylines or splines)
 250     isConnectedBySegment: boolean, to add a segment connecting the polylines, default False
 251     tolerance: used by the algorithm to detect coincidences, default 1.E-3
 252     return:
 253     the merged polyline
 254     """
 255     op=HYDROData_PolylineOperator()
 256     op.Merge(document, polyName, polyLines, isConnectedBySegment, tolerance)
 257     shape = document.FindObjectByName(polyName)
 258     return shape
 259
 260 def createAxis3DDEmbankmentAltiProfile(document, axis3DName, axisXY, altiPts):
 261     """
 262     Create a 3D polyline to be used as an axis for an embankement, using a 2D axis and a list of altitudes along the axis.
 263     parameters:
 264     document: current HYDROData_Document
 265     axis3DName: name to give to the axis3D created
 266     axisXY: 2D polyline (spline), Embankment line
 267     altiPts: list of tuples(x, h) : parametric coordinate, altitude, along the axisXY
 268     return:
 269     altiProfile, axis3D
 270     """
 271     altiProfile = document.CreateObject(KIND_PROFILE)
 272     altiProfile.SetName("%s_altiProfile"%axis3DName)
 273     altiProfilePoints = []
 274     for p in altiPts:
 275         altiProfilePoints.append(gp_XY(p[0], p[1]))
 276     altiProfile.SetParametricPoints(altiProfilePoints)
 277     altiProfile.Update()
 278     axis3D = document.CreateObject(KIND_POLYLINE)
 279     axis3D.SetName(axis3DName)
 280     axis3D.SetPolylineXY(axisXY)
 281     axis3D.SetProfileUZ(altiProfile.GetProfileUZ())
 282     axis3D.Update()
 283     return (altiProfile, axis3D)
 284
 285 def createAxis3DDEmbankmentBathy(document, axis3DName, axisXY, aCloud):
 286     """
 287     Create a 3D polyline to be used as an axis for an embankement, using a 2D axis and a bathymetry to project the axis.
 288     parameters:
 289     document: current HYDROData_Document
 290     axis3DName: name to give to the axis3D created
 291     axisXY: imported polyline, Embankment line
 292     aCloud:  altimetry, for projection of axisXY
 293     return:
 294     altiProfile, axis3D
 295     """
 296     altiProfile = document.CreateObject(KIND_PROFILE)
 297     altiProfile.SetName("%s_altiProfile"%axis3DName)
 298     axis3D = document.CreateObject(KIND_POLYLINE)
 299     axis3D.SetName(axis3DName)
 300     axis3D.SetPolylineXY(axisXY)
 301     axis3D.SetAltitudeObject(aCloud)
 302     axis3D.SetChildProfileUZ(altiProfile.GetProfileUZ())
 303     altiProfile.Update()
 304     axis3D.Update()
 305     return (altiProfile, axis3D)
 306
 307 def createEmbankmentSectionA(document, embankmentName, axis3D, sectionPoints, d, displayLevel):
 308     """
 309     Define the section of an embankement and extrude it along a 3D axis. Section is defined with a list of points (axis distance, altitude).
 310     parameters:
 311     document: current HYDROData_Document
 312     embankmentName: name to give to the embankment
 313     axis3D: the 3D axis of the embankment
 314     sectionPoints: a list of tuple (x, h) to define a half section, beginning to x=0 (center)
 315     d calculation parameter, take 2 or 3 times the section width
 316     displayLevel : z level for 2D representation: high values are drawn above low values
 317     return:
 318     section, embankment
 319     """
 320     section = document.CreateObject(KIND_PROFILE)
 321     section.SetName("%s_section"%embankmentName)
 322     sectionPts = []
 323     for i,p in enumerate(sectionPoints):
 324         sectionPts.append(gp_XY(p[0], p[1]))
 325         if i>0:
 326             sectionPts.insert(0, gp_XY(-p[0], p[1]))
 327     section.SetParametricPoints(sectionPts)
 328     section.Update()
 329     embankment = document.CreateObject(KIND_DIGUE)
 330     embankment.SetName(embankmentName)
 331     embankment.SetGuideLine(axis3D)
 332     embankment.SetProfileMode(True)
 333     embankment.SetProfile(section)
 334     embankment.SetEquiDistance(d)
 335     embankment.SetZLevel(displayLevel)
 336     embankment.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 337     embankment.Update()
 338     return (section, embankment)
 339
 340 def createEmbankmentSectionB(document, embankmentName, axis3D, LC,DZ,CZ, d, displayLevel):
 341     """
 342     Define the section of an embankement and extrude it along a 3D axis. Section is defined with a width, an height, the position of the axis.
 343     parameters:
 344     document: current HYDROData_Document
 345     embankmentName: name to give to the embankment
 346     axis3D: the 3D axis of the embankment
 347     LC, DZ, CZ: width, height, height offset above axis3D
 348     d calculation parameter, take 2 or 3 times the section width
 349     displayLevel : z level for 2D representation: high values are drawn above low values
 350     return:
 351     embankment
 352     """
 353     embankment = document.CreateObject(KIND_DIGUE)
 354     embankment.SetName(embankmentName)
 355     embankment.SetGuideLine(axis3D)
 356     embankment.SetProfileMode(False)
 357     embankment.SetLCValue(LC)
 358     embankment.SetDeltaZValue(DZ)
 359     embankment.SetCoteZValue(CZ)
 360     embankment.SetEquiDistance(d)
 361     embankment.SetZLevel(displayLevel)
 362     embankment.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 363     embankment.Update()
 364     return embankment
 365
 366 def loadStricklerTable(document, stricklerFile):
 367     """
 368     Load a table of Corine Land Cover types  and codes, used for Strickler coefficients (editable text format)
 369     parameters:
 370     document: current HYDROData document
 371     stricklerFile: full path of the table (*.txt)
 372     return:
 373     table loaded, code corresponding to Corine Land Cover used for the table ('CODE_06', 'CODE_12'...)
 374     """
 375     code = ""
 376     with open(stricklerFile) as f:
 377         line = f.readline()
 378         code = line.split()[0]
 379     a = os.path.splitext(stricklerFile)
 380     stricklerName = os.path.basename(a[0])
 381     strickler_table = document.CreateObject( KIND_STRICKLER_TABLE )
 382     strickler_table.Import(stricklerFile)
 383     strickler_table.SetName(stricklerName)
 384     strickler_table.SetAttrName(code)
 385     strickler_table.Update()
 386     return strickler_table, code
 387
 388 def loadLandCoverMap(document, landCoverFile, stricklerFile, displayLevel):
 389     """
 390     Load a Corine Land Cover Map
 391     parameters:
 392     document: current HYDROData document
 393     landCoverFile: shapeFile defining the map (*.shp)
 394     stricklerFile: full path of the table (*.txt)
 395     displayLevel : z level for 2D representation: high values are drawn above low values
 396     return:
 397     land cover map loaded
 398     """
 399     alltypes = []
 400     allvalues = []
 401     code = ""
 402     l = 0
 403     with open(stricklerFile) as f:
 404         for line in f:
 405             items = line.split('"')
 406             if l ==0:
 407                 code = line.split()[0]
 408             else:
 409                 alltypes.append(items[1])
 410                 itm2 = items[-1].split()
 411                 allvalues.append(itm2[-1])
 412             l = l + 1
 413     print(alltypes)
 414     print(allvalues)
 415     while code[0] != 'C':    # the code begins with 'C'
 416         code = code[1:]      # get rid of encoding
 417     print (code)
 418
 419     a = os.path.splitext(landCoverFile)
 420     landCoverName = os.path.basename(a[0])
 421     landCoverDB = a[0] + '.dbf'
 422     landCover = document.CreateObject(KIND_LAND_COVER_MAP)
 423     print(landCoverFile)
 424     print(landCoverDB)
 425     print(landCoverName)
 426     landCover.SetName(landCoverName)
 427     landCover.SetZLevel(displayLevel)
 428     if not landCover.ImportSHP(landCoverFile):
 429         raise ValueError('problem while loading LandCoverMap shape')
 430     if landCover.ImportDBF(landCoverDB, code, allvalues, alltypes ) != landCover.DBFStatus_OK:
 431         raise ValueError('problem while loading LandCoverMap data base')
 432     landCover.Update()
 433     return landCover
 434
 435
 436
 437
 438
 439