From d32791da491a27db3766ada199d9066b51b1dd7e Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Wed, 12 Jul 2017 16:40:20 +0200 Subject: [PATCH] corrections sur les scripts avec MedCoupling --- src/HYDROTools/interpolS.py | 176 +++++++++++++++++++++--------------- src/HYDROTools/interpolZ.py | 61 +++++++++---- 2 files changed, 145 insertions(+), 92 deletions(-) mode change 100644 => 100755 src/HYDROTools/interpolS.py diff --git a/src/HYDROTools/interpolS.py b/src/HYDROTools/interpolS.py old mode 100644 new mode 100755 index 3376db44..55c017a6 --- a/src/HYDROTools/interpolS.py +++ b/src/HYDROTools/interpolS.py @@ -1,90 +1,120 @@ # -*- coding: utf-8 -*- +""" +assignStrickler +""" +__revision__ = "V2.02" + import os -import sys -import shutil import salome -from med import medfile -from med import medmesh -from med import medfield -from med import medenum - import HYDROPy +import numpy as np +import MEDLoader as ml +import MEDCoupling as mc + # Get current study id salome.salome_init() theStudy = salome.myStudy theStudyId = salome.myStudyId -"""Bad parameters exception""" class BadParamsError(ValueError): + """Bad parameters exception""" pass -""" -Creates MED file with scalar field associated with the mesh of the calculation case. -The scalar value associated with the mesh node is the Strickler coefficient of the land cover containing the node. -case_name: calculation case name in the study -med_file_name: path to input MED file with mesh on nodes -output_file_name: path to output MED file with -med_field_name: field name -""" -def assignStrickler(case_name, med_file_name, output_file_name, med_field_name='BOTTOM FRICTION'): - # Check calculation case - doc = HYDROPy.HYDROData_Document.Document( theStudyId ) - case = doc.FindObjectByName(case_name) - if case is None: - raise BadParamsError("Calculation case '%s' not found" % case_name) - - # Check input MED file - if not os.path.exists(med_file_name): - raise BadParamsError("Input file '%s' not exists" % med_file_name) - - # Copy input file to the output path, if the output path is empty - use the input path for the output - file_path = med_file_name - if output_file_name and (output_file_name != file_path): - shutil.copyfile(med_file_name, output_file_name) - file_path = output_file_name - - # Open input MED file - fid = medfile.MEDfileOpen(file_path, medenum.MED_ACC_RDEXT) - - # Get mesh info - mesh_name, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1) - - nb_nodes, chgt, trsf = medmesh.MEDmeshnEntity(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT, - medenum.MED_NODE, medenum.MED_NONE, - medenum.MED_COORDINATE, medenum.MED_NO_CMODE) - - - # Create field - if nb_nodes > 0: - # Get node coordinates - coords = medfile.MEDFLOAT(nb_nodes*sdim) - medmesh.MEDmeshNodeCoordinateRd(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT, - medenum.MED_FULL_INTERLACE, coords) - - # Get list of Strickler coefficient values for the nodes - values = [] - x_coords = coords[0::sdim] - y_coords = coords[1::sdim] - #print "case.GetStricklerCoefficientForPoints" + +def assignStrickler(nomCas, input_file_name, output_file_name="", med_field_name='BOTTOM FRICTION', verbose=False): + """ + assignStrickler creates the scalar field associated with the mesh of the calculation case that + represents the Strickler coefficient of the land cover. + + In: + nomCas: Calculation Case Name in module HYDRO + input_file_name: med file name corresponding to the HYDRO case + output_file_name: med file name with the Strickler coefficient; default is the same as fichierMaillage + med_field_name: name of the field of the Strickler coefficient + default value is 'BOTTOM FRICTION'. + Out: + """ + erreur = 0 + message = "" + + while not erreur: + + if verbose: + print "nomCas:", nomCas + print "input_file_name:", input_file_name + print "output_file_name:", output_file_name + print "med_field_name:", med_field_name + +# 1. Controls +# 1.1. Check calculation case + doc = HYDROPy.HYDROData_Document.Document(theStudyId) + case = doc.FindObjectByName(nomCas) + if case is None: + raise BadParamsError("Calculation case '%s' not found" % nomCas) + +# 1.2. Check input MED file + if not os.path.exists(input_file_name): + raise BadParamsError("Input file '%s' not exists" % input_file_name) + +# 2. The output file +# 2.1. If output MED file is not defined, it is equal to the input file + if not output_file_name: + output_file_name = input_file_name + +# 2.2. Copy input file to the output path, if the output path is different from the input path + if output_file_name != input_file_name: + import shutil + shutil.copyfile(input_file_name, output_file_name) +# +# 3. Reads the mesh +# + meshMEDFileRead = ml.MEDFileMesh.New(output_file_name) +# +# 4. Gets the information about the nodes +# + nbnodes = meshMEDFileRead.getNumberOfNodes() + if verbose: + print "Number of nodes: %d" % nbnodes +# + coords = meshMEDFileRead.getCoords() + #print "coords =\n", coords +# +# 5. Calculation of the Strickler coefficient + values = list() + x_coords = coords[:, 0] + y_coords = coords[:, 1] + #print "x_coords =\n", x_coords + #print "y_coords =\n", y_coords values = case.GetStricklerCoefficientForPoints(x_coords, y_coords, 0.0, True) - - # Write the values to the field - if len(values) > 0: - values = medfile.MEDFLOAT(values) - - comp = "strickler coeff " - medfield.MEDfieldCr(fid, med_field_name, medfile.MED_FLOAT64, - 1, comp, '', '', mesh_name) - - medfield.MEDfieldValueWr(fid, med_field_name, 0, medenum.MED_NO_IT, 0.0, - medenum.MED_NODE, medenum.MED_NONE, - medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, len(values), values) - - print "MED field '%s' on %s nodes has been created." % (med_field_name, nb_nodes) - - # close MED file - medfile.MEDfileClose(fid) +# +# 6. Field MED file +# + coeff = mc.DataArrayDouble(np.asfarray(values, dtype='float')) + coeff.setInfoOnComponents(["Strickler [SI]"]) + #print "coeff =\n", coeff + if verbose: + print ".. Ecriture du Strickler sous le nom '"+med_field_name+"'" + fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES) + fieldOnNodes.setName(med_field_name) + fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0)) + fieldOnNodes.setArray(coeff) +# Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC +# instant = 0.0 +# numero d'itération : 0 +# pas de numero d'ordre (-1) + fieldOnNodes.setTime(0.0, 0, -1) +# + fMEDFile_ch_d = ml.MEDFileField1TS() + fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes) + fMEDFile_ch_d.write(output_file_name, 0) +# + break + + if erreur: + print message + + return diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py index 94618d72..318ec829 100755 --- a/src/HYDROTools/interpolZ.py +++ b/src/HYDROTools/interpolZ.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """ -example of use case of interpolZ and createZfield methods: +example of use case of interpolZ: # --- case name in HYDRO nomCas = 'inondation1' @@ -20,7 +20,7 @@ zUndef = 90 # --- Z interpolation on the bathymety/altimetry on the mesh nodes interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod) """ -__revision__ = "V2.02" +__revision__ = "V2.04" # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -46,6 +46,7 @@ from med import medenum #from med import medlocalization #from med import medlink +# La fonction createZfield1 ne sert plus à partir du 12/07/2017 def createZfield1(fichierMaillage): """ Complete the mesh for Telemac. @@ -108,9 +109,7 @@ def createZfield1(fichierMaillage): # ----------------------------------------------------------------------------- -import MEDLoader -from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh - +# La fonction createZfield2 ne sert plus à partir du 12/07/2017 def createZfield2(fichierMaillage): """ Complete the mesh for Telemac. @@ -121,6 +120,9 @@ def createZfield2(fichierMaillage): return F.med : med file containing the field "BOTTOM" """ + import MEDLoader + from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh + basename = fichierMaillage[:-4] fichierZMaillage = basename + 'Z.med' fichierFMaillage = basename + 'F.med' @@ -152,7 +154,7 @@ class MyInterpolator( HYDROPy.HYDROData_IInterpolator ): # ----------------------------------------------------------------------------- -def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, zUndefThreshold=-9000., xyzFile=False, verbose=False): +def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False): """ interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node. @@ -161,19 +163,19 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet nomCas: Calculation Case Name in module HYDRO fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions. - Key: face group name, value: region name in the HYDRO Case + Key: face group name + Value: region name in the HYDRO Case zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct). default value is 90. interpolMethod: integer value 0 = nearest point on bathymetry (default) 1 = linear interpolation zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct). - zUndefThreshold: if the interpolated z is under this value, it is replaced by zUndef - default value is -9000. xyzFile: True/False to write an ascii file with xyz for every node. Default is False. Out: statz: statistique for z - Key: face group name, value: (minz, maxz, meanz, stdz, v05z, v95z) + Key: face group name + Value: (minz, maxz, meanz, stdz, v05z, v95z) Out: return F.med : med file with Z value on nodes and in a field "BOTTOM" """ @@ -182,11 +184,11 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet message = "" while not erreur: - + if verbose: + print "nomCas:", nomCas print "interpolMethod: %d" % interpolMethod print "zUndef:", zUndef - print "zUndefThreshold: %f" % zUndefThreshold doc = HYDROPy.HYDROData_Document.Document(theStudyId) cas = doc.FindObjectByName(nomCas) @@ -203,7 +205,7 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet fichierFonds = basename + '.xyz' print "fichierFonds =", fichierFonds # -# 1. Reads the mesh and gets the number of nodes +# 1. Reads the mesh # meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage) # @@ -315,7 +317,17 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet bathy[nodeId] = vz[iaux] # # 5. Minimum: -# If the value is lower than a threshold, an "undefined" valeur is set +# During the interpolation, if no value is available over a node, a default value +# is set: -9999. It has no importance for the final computation, but if the field +# or the mesh is displayed, it makes huge gap. To prevent this artefact, a more +# convenient "undefined" value is set. This new undefinde value is given by the user. +# +# zUndefThreshold: the default value is -9000. It is tied with the value -9999. given +# by the interpolation when no value is defined. +# + zUndefThreshold = -9000. + if verbose: + print "zUndefThreshold: %f" % zUndefThreshold# # #print "bathy :\n", bathy np_aux_z = (bathy < zUndefThreshold).nonzero() @@ -341,8 +353,12 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet # 7.1. Modification of the z coordinates # bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float')) - coords3D = DataArrayDouble.Meld([coords, bathy_dd]) - #print "coords3D\n", coords3D + bathy_dd.setInfoOnComponents(["Z [m]"]) +# + coords3D = ml.DataArrayDouble.Meld([coords, bathy_dd]) + coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"]) + #print "coords3D =\n", coords3D +# meshMEDFileRead.setCoords(coords3D) # # 7.2. Writes the 3D mesh @@ -353,12 +369,19 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet # # 7.3. Writes the field # + med_field_name = "BOTTOM" if verbose: - print ".. Ecriture du champ BOTTOM" - fieldOnNodes=ml.MEDCouplingFieldDouble(ml.ON_NODES, ml.ONE_TIME) - fieldOnNodes.setName("BOTTOM") + print ".. Ecriture du champ '"+med_field_name+"'" + #print "bathy_dd =\n", bathy_dd + fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES) + fieldOnNodes.setName(med_field_name) fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0)) fieldOnNodes.setArray(bathy_dd) +# Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC +# instant = 0.0 +# numero d'itération : 0 +# pas de numero d'ordre (-1) + fieldOnNodes.setTime(0.0, 0, -1) # fMEDFile_ch_d = ml.MEDFileField1TS() fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes) -- 2.39.2