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, shapeName, 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     shapeName: used to retrieve the objects in the HYDRO document: see above, name is given without the sequence suffix
  91     isSpline: boolean to convert the polylines in  splines
  92     displayLevel: integer >=0, the higher levels are displayed on top
  93     return:
  94     list of shape objects
  95     """
  96     HYDROData_PolylineXY.ImportShapesFromFile(shapeFile)
  97     shapeType = 0 # polyline
  98     if iSpline:
  99         shapeType = 1 # polyline
 100     shapes = []
 101     isShapeFound = True
 102     print("importPolylines %s %s"%(shapeName, shapeType))
 103     index = 0  # for multiple shapes
 104     while isShapeFound:
 105         shapeNameIndex = shapeName + "_%d" % index
 106         index = index + 1
 107         print("try name %s"%shapeNameIndex)
 108         shape = document.FindObjectByName(shapeNameIndex)
 109         if shape is None:
 110             isShapeFound = False
 111         else:
 112             print("found %s"%shapeNameIndex)
 113             for i in range(shape.NbSections()):
 114                 shape.SetSectionType(i, shapeType)
 115                 shape.Update()
 116                 shape.SetZLevel(displayLevel)
 117             shapes.append(shape)
 118     if index <= 1: # no multiple shapes, try single shape
 119         print("try name %s"%shapeName)
 120         shape = document.FindObjectByName(shapeName) # single shape
 121         if shape is not None:
 122             print("found %s"%shapeName)
 123             for i in range(shape.NbSections()):
 124                 shape.SetSectionType(i, shapeType)
 125                 shape.Update()
 126                 shape.SetZLevel(displayLevel)
 127             shapes.append(shape)
 128     return shapes
 129
 130 def splitShapeTool(document, aShape, toolShape, tolerance=1.e-2):
 131     """
 132     Split the shape (2D polyline or spline) by the shape tool.
 133     return a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 134     parameters;
 135     document: current HYDROData document
 136     aShape: a loaded shape (2D polyline or spline)
 137     toolShape: a loaded shape (2D polyline or spline)
 138     tolerance: used by the algorithm to detect intersections, default value 1.e-2
 139     return:
 140     a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 141     """
 142     name = aShape.GetName()
 143     existingNbSplit = {}
 144     index = 1
 145     found = True
 146     while found:
 147         nameIndex = name + "_%d"%(index)
 148         shape = document.FindObjectByName(nameIndex)
 149         if shape is None:
 150             found = False
 151         else:
 152             index = index + 1
 153     existingNbSplit[name] = index
 154     print(existingNbSplit)
 155
 156     op = HYDROData_PolylineOperator()
 157     op.SplitTool(document, aShape, toolShape, tolerance)
 158
 159     listSplit = []
 160     index = existingNbSplit[name]
 161     found = True
 162     while found:
 163         nameIndex = name + "_%d"%(index)
 164         shape = document.FindObjectByName(nameIndex)
 165         if shape is None:
 166             found = False
 167         else:
 168             print("found %s"%nameIndex)
 169             index = index + 1
 170             listSplit.append(shape)
 171     return listSplit
 172
 173 def splitShapesAll(document, shapeslist, precision=1.E-2):
 174     """
 175     Split all the shapes (2D polylines or splines) in the list by all the other shapes.
 176     return a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 177     parameters;
 178     document: current HYDROData document
 179     shapeslist: a list of loaded shapes (2D polylines or splines)
 180     precision: used by the algorithm to detect intersections, default value 1.e-2
 181     return:
 182     a list of all the split shapes, named after their original shape name, with a sequence suffix _n
 183     """
 184     names = [s.GetName() for s in shapeslist]
 185
 186     existingNbSplit = {}
 187     for name in names:
 188         index = 1
 189         found = True
 190         while found:
 191             nameIndex = name + "_%d"%(index)
 192             shape = document.FindObjectByName(nameIndex)
 193             if shape is None:
 194                 found = False
 195             else:
 196                 index = index + 1
 197         existingNbSplit[name] = index
 198     print(existingNbSplit)
 199
 200     op = HYDROData_PolylineOperator()
 201     op.SplitAll(document, shapeslist, precision)
 202
 203     listSplit = []
 204     for name in names:
 205         index = existingNbSplit[name]
 206         found = True
 207         while found:
 208             nameIndex = name + "_%d"%(index)
 209             shape = document.FindObjectByName(nameIndex)
 210             if shape is None:
 211                 found = False
 212             else:
 213                 print("found %s"%nameIndex)
 214                 index = index + 1
 215                 listSplit.append(shape)
 216     return listSplit
 217
 218
 219 def importBathymetry(document, bathyFile):
 220     """
 221     import a bathymetry file (*.xyz or *.asc)
 222     parameters:
 223     document: current HYDROData document
 224     bathyFile: full path of the bathymetry (*.xyz or *.shp)
 225     return:
 226     the bathymetry object
 227     """
 228     bathy = document.CreateObject(KIND_BATHYMETRY)
 229     a = os.path.splitext(bathyFile)
 230     bathyName = os.path.basename(a[0])
 231     bathy.SetName(bathyName)
 232     bathy.SetAltitudesInverted(0)
 233     if not(bathy.ImportFromFile(bathyFile)):
 234         raise ValueError('problem while loading bathymetry')
 235     bathy.Update()
 236     return bathy
 237
 238 def createImmersibleZone(document, imzName, polyLine, bathy, isImmersible, displayLevel):
 239     """
 240     Create an immersible or unsubmersible zone with a closed polyline and a bathymetry.
 241     parameters:
 242     document: current HYDROData document
 243     imzName: the name to give to the zone object
 244     polyline: the closed 2D polyline or spline delimiting the zone
 245     bathy: the bathymetry to associate to the zone
 246     isImmersible: boolean
 247     displayLevel: integer >=0, the higher levels are displayed on top
 248     return:
 249     the zone created
 250     """
 251     imz = document.CreateObject(KIND_IMMERSIBLE_ZONE)
 252     imz.SetName(imzName)
 253     imz.SetZLevel(displayLevel)
 254     if bathy is not None:
 255         imz.SetAltitudeObject(bathy)
 256     imz.SetPolyline(polyLine)
 257     imz.SetIsSubmersible(isImmersible)
 258     imz.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 259     imz.Update()
 260     return imz
 261
 262 def mergePolylines(document, polyName, polyLines, isConnectedBySegment = False, tolerance = 1.E-3):
 263     """
 264     Regroup several 2D polylines or spline in one, with several sections
 265     parameters:
 266     document: current HYDROData document
 267     polyName: name to give to the new object
 268     polyLines: a list of loaded shapes (2D polylines or splines)
 269     isConnectedBySegment: boolean, to add a segment connecting the polylines, default False
 270     tolerance: used by the algorithm to detect coincidences, default 1.E-3
 271     return:
 272     the merged polyline
 273     """
 274     op=HYDROData_PolylineOperator()
 275     op.Merge(document, polyName, polyLines, isConnectedBySegment, tolerance)
 276     shape = document.FindObjectByName(polyName)
 277     return shape
 278
 279 def createAxis3DDEmbankmentAltiProfile(document, axis3DName, axisXY, altiPts):
 280     """
 281     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.
 282     parameters:
 283     document: current HYDROData_Document
 284     axis3DName: name to give to the axis3D created
 285     axisXY: 2D polyline (spline), Embankment line
 286     altiPts: list of tuples(x, h) : parametric coordinate, altitude, along the axisXY
 287     return:
 288     altiProfile, axis3D
 289     """
 290     altiProfile = document.CreateObject(KIND_PROFILE)
 291     altiProfile.SetName("%s_altiProfile"%axis3DName)
 292     altiProfilePoints = []
 293     for p in altiPts:
 294         altiProfilePoints.append(gp_XY(p[0], p[1]))
 295     altiProfile.SetParametricPoints(altiProfilePoints)
 296     altiProfile.Update()
 297     axis3D = document.CreateObject(KIND_POLYLINE)
 298     axis3D.SetName(axis3DName)
 299     axis3D.SetPolylineXY(axisXY)
 300     axis3D.SetProfileUZ(altiProfile.GetProfileUZ())
 301     axis3D.Update()
 302     return (altiProfile, axis3D)
 303
 304 def createAxis3DDEmbankmentBathy(document, axis3DName, axisXY, aCloud):
 305     """
 306     Create a 3D polyline to be used as an axis for an embankement, using a 2D axis and a bathymetry to project the axis.
 307     parameters:
 308     document: current HYDROData_Document
 309     axis3DName: name to give to the axis3D created
 310     axisXY: imported polyline, Embankment line
 311     aCloud:  altimetry, for projection of axisXY
 312     return:
 313     altiProfile, axis3D
 314     """
 315     altiProfile = document.CreateObject(KIND_PROFILE)
 316     altiProfile.SetName("%s_altiProfile"%axis3DName)
 317     axis3D = document.CreateObject(KIND_POLYLINE)
 318     axis3D.SetName(axis3DName)
 319     axis3D.SetPolylineXY(axisXY)
 320     axis3D.SetAltitudeObject(aCloud)
 321     axis3D.SetChildProfileUZ(altiProfile.GetProfileUZ())
 322     altiProfile.Update()
 323     axis3D.Update()
 324     return (altiProfile, axis3D)
 325
 326 def createEmbankmentSectionA(document, embankmentName, axis3D, sectionPoints, d, displayLevel):
 327     """
 328     Define the section of an embankement and extrude it along a 3D axis. Section is defined with a list of points (axis distance, altitude).
 329     parameters:
 330     document: current HYDROData_Document
 331     embankmentName: name to give to the embankment
 332     axis3D: the 3D axis of the embankment
 333     sectionPoints: a list of tuple (x, h) to define a half section, beginning to x=0 (center)
 334     d calculation parameter, take 2 or 3 times the section width
 335     displayLevel : z level for 2D representation: high values are drawn above low values
 336     return:
 337     section, embankment
 338     """
 339     section = document.CreateObject(KIND_PROFILE)
 340     section.SetName("%s_section"%embankmentName)
 341     sectionPts = []
 342     for i,p in enumerate(sectionPoints):
 343         sectionPts.append(gp_XY(p[0], p[1]))
 344         if i>0:
 345             sectionPts.insert(0, gp_XY(-p[0], p[1]))
 346     section.SetParametricPoints(sectionPts)
 347     section.Update()
 348     embankment = document.CreateObject(KIND_DIGUE)
 349     embankment.SetName(embankmentName)
 350     embankment.SetGuideLine(axis3D)
 351     embankment.SetProfileMode(True)
 352     embankment.SetProfile(section)
 353     embankment.SetEquiDistance(d)
 354     embankment.SetZLevel(displayLevel)
 355     embankment.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 356     embankment.Update()
 357     return (section, embankment)
 358
 359 def createEmbankmentSectionB(document, embankmentName, axis3D, LC,DZ,CZ, d, displayLevel):
 360     """
 361     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.
 362     parameters:
 363     document: current HYDROData_Document
 364     embankmentName: name to give to the embankment
 365     axis3D: the 3D axis of the embankment
 366     LC, DZ, CZ: width, height, height offset above axis3D
 367     d calculation parameter, take 2 or 3 times the section width
 368     displayLevel : z level for 2D representation: high values are drawn above low values
 369     return:
 370     embankment
 371     """
 372     embankment = document.CreateObject(KIND_DIGUE)
 373     embankment.SetName(embankmentName)
 374     embankment.SetGuideLine(axis3D)
 375     embankment.SetProfileMode(False)
 376     embankment.SetLCValue(LC)
 377     embankment.SetDeltaZValue(DZ)
 378     embankment.SetCoteZValue(CZ)
 379     embankment.SetEquiDistance(d)
 380     embankment.SetZLevel(displayLevel)
 381     embankment.SetFillingColor(QColor(randint(0,255), randint(0,255), randint(0,255), 255))
 382     embankment.Update()
 383     return embankment
 384
 385 def loadStricklerTable(document, stricklerFile):
 386     """
 387     Load a table of Corine Land Cover types  and codes, used for Strickler coefficients (editable text format)
 388     parameters:
 389     document: current HYDROData document
 390     stricklerFile: full path of the table (*.txt)
 391     return:
 392     table loaded, code corresponding to Corine Land Cover used for the table ('CODE_06', 'CODE_12'...)
 393     """
 394     code = ""
 395     with open(stricklerFile) as f:
 396         line = f.readline()
 397         code = line.split()[0]
 398     a = os.path.splitext(stricklerFile)
 399     stricklerName = os.path.basename(a[0])
 400     strickler_table = document.CreateObject( KIND_STRICKLER_TABLE )
 401     strickler_table.Import(stricklerFile)
 402     strickler_table.SetName(stricklerName)
 403     strickler_table.SetAttrName(code)
 404     strickler_table.Update()
 405     return strickler_table, code
 406
 407 def loadLandCoverMap(document, landCoverFile, stricklerFile, displayLevel):
 408     """
 409     Load a Corine Land Cover Map
 410     parameters:
 411     document: current HYDROData document
 412     landCoverFile: shapeFile defining the map (*.shp)
 413     stricklerFile: full path of the table (*.txt)
 414     displayLevel : z level for 2D representation: high values are drawn above low values
 415     return:
 416     land cover map loaded
 417     """
 418     alltypes = []
 419     allvalues = []
 420     code = ""
 421     l = 0
 422     with open(stricklerFile) as f:
 423         for line in f:
 424             items = line.split('"')
 425             if l ==0:
 426                 code = line.split()[0]
 427             else:
 428                 alltypes.append(items[1])
 429                 itm2 = items[-1].split()
 430                 allvalues.append(itm2[-1])
 431             l = l + 1
 432     print(alltypes)
 433     print(allvalues)
 434     while code[0] != 'C':    # the code begins with 'C'
 435         code = code[1:]      # get rid of encoding
 436     print (code)
 437
 438     a = os.path.splitext(landCoverFile)
 439     landCoverName = os.path.basename(a[0])
 440     landCoverDB = a[0] + '.dbf'
 441     landCover = document.CreateObject(KIND_LAND_COVER_MAP)
 442     print(landCoverFile)
 443     print(landCoverDB)
 444     print(landCoverName)
 445     landCover.SetName(landCoverName)
 446     landCover.SetZLevel(displayLevel)
 447     if not landCover.ImportSHP(landCoverFile):
 448         raise ValueError('problem while loading LandCoverMap shape')
 449     if landCover.ImportDBF(landCoverDB, code, allvalues, alltypes ) != landCover.DBFStatus_OK:
 450         raise ValueError('problem while loading LandCoverMap data base')
 451     landCover.Update()
 452     return landCover
 453
 454
 455
 456
 457
 458