From 2650061db024abb253f34cdc26685433d6e426d4 Mon Sep 17 00:00:00 2001 From: mzn Date: Fri, 16 Sep 2016 15:33:47 +0300 Subject: [PATCH] Lot5: Add script for creation of field with Strickler coefficient values. --- src/HYDROData/HYDROData_CalculationCase.cxx | 6 ++ src/HYDROPy/CAS/vector.sip | 16 ++-- src/HYDROPy/HYDROData_CalculationCase.sip | 31 +++++++ src/HYDROTools/CMakeLists.txt | 1 + src/HYDROTools/interpolS.py | 89 +++++++++++++++++++++ 5 files changed, 135 insertions(+), 8 deletions(-) create mode 100644 src/HYDROTools/interpolS.py diff --git a/src/HYDROData/HYDROData_CalculationCase.cxx b/src/HYDROData/HYDROData_CalculationCase.cxx index fef44bd1..f82e4e32 100644 --- a/src/HYDROData/HYDROData_CalculationCase.cxx +++ b/src/HYDROData/HYDROData_CalculationCase.cxx @@ -171,6 +171,12 @@ QStringList HYDROData_CalculationCase::DumpToPython( const QString& thePyS Handle(HYDROData_PolylineXY) aBoundaryPolyline = GetBoundaryPolyline(); setPythonReferenceObject( thePyScriptPath, theTreatedObjects, aResList, aBoundaryPolyline, "SetBoundaryPolyline" ); + Handle(HYDROData_StricklerTable) aStricklerTable = GetStricklerTable(); + setPythonReferenceObject( thePyScriptPath, theTreatedObjects, aResList, aStricklerTable, "SetStricklerTable" ); + + Handle(HYDROData_LandCoverMap) aLandCoverMap = GetLandCoverMap(); + setPythonReferenceObject( thePyScriptPath, theTreatedObjects, aResList, aLandCoverMap, "SetLandCoverMap" ); + if( aMode==AUTOMATIC ) DumpRulesToPython( aCalculName, aResList ); diff --git a/src/HYDROPy/CAS/vector.sip b/src/HYDROPy/CAS/vector.sip index 8ea02b9c..d3823364 100644 --- a/src/HYDROPy/CAS/vector.sip +++ b/src/HYDROPy/CAS/vector.sip @@ -36,7 +36,7 @@ template return NULL; // Set the list elements. - for ( int i = 1, n = sipCpp->size(); i <= n; ++i ) + for ( int i = 0, n = sipCpp->size(); i < n; ++i ) { TYPE* t = new TYPE( (*sipCpp)[i] ); @@ -49,7 +49,7 @@ template return NULL; } - PyList_SET_ITEM( l, i - 1, pobj ); + PyList_SET_ITEM( l, i, pobj ); } return l; @@ -124,7 +124,7 @@ template return NULL; // Set the list elements. - for ( int i = 1, n = sipCpp->size(); i <= n; ++i ) + for ( int i = 0, n = sipCpp->size(); i < n; ++i ) { PyObject *pobj; if ( ( pobj = PyFloat_FromDouble( (*sipCpp)[i] ) ) == NULL ) @@ -134,7 +134,7 @@ template return NULL; } - PyList_SET_ITEM( l, i - 1, pobj ); + PyList_SET_ITEM( l, i, pobj ); } return l; @@ -186,7 +186,7 @@ template return NULL; // Set the list elements. - for ( int i = 1, n = sipCpp->size(); i <= n; ++i ) + for ( int i = 0, n = sipCpp->size(); i < n; ++i ) { PyObject *pobj; if ( ( pobj = SIPLong_FromLong( (*sipCpp)[i] ) ) == NULL ) @@ -196,7 +196,7 @@ template return NULL; } - PyList_SET_ITEM( l, i - 1, pobj ); + PyList_SET_ITEM( l, i, pobj ); } return l; @@ -248,7 +248,7 @@ template return NULL; // Set the list elements. - for ( int i = 1, n = sipCpp->size(); i <= n; ++i ) + for ( int i = 0, n = sipCpp->size(); i < n; ++i ) { PyObject *pobj; if ( ( pobj = SIPLong_FromLong( (*sipCpp)[i] ) ) == NULL ) @@ -258,7 +258,7 @@ template return NULL; } - PyList_SET_ITEM( l, i - 1, pobj ); + PyList_SET_ITEM( l, i, pobj ); } return l; diff --git a/src/HYDROPy/HYDROData_CalculationCase.sip b/src/HYDROPy/HYDROData_CalculationCase.sip index 01b1ae04..1bbc51e6 100644 --- a/src/HYDROPy/HYDROData_CalculationCase.sip +++ b/src/HYDROPy/HYDROData_CalculationCase.sip @@ -215,7 +215,38 @@ public: * Remove reference boundary polyline object from calculation case. */ void RemoveBoundaryPolyline(); + + void SetLandCoverMap( HYDROData_LandCoverMap theLandCoverMap ) [void ( const Handle_HYDROData_LandCoverMap& )]; + %MethodCode + Handle(HYDROData_LandCoverMap) aRef = + Handle(HYDROData_LandCoverMap)::DownCast( createHandle( a0 ) ); + if ( !aRef.IsNull() ) + { + Py_BEGIN_ALLOW_THREADS + if ( sipSelfWasArg ) { + sipCpp->HYDROData_CalculationCase::SetLandCoverMap( aRef ); + } else { + sipCpp->SetLandCoverMap( aRef ); + } + Py_END_ALLOW_THREADS + } + %End + void SetStricklerTable( HYDROData_StricklerTable theStricklerTable ) [void ( const Handle_HYDROData_StricklerTable& )]; + %MethodCode + Handle(HYDROData_StricklerTable) aRef = + Handle(HYDROData_StricklerTable)::DownCast( createHandle( a0 ) ); + if ( !aRef.IsNull() ) + { + Py_BEGIN_ALLOW_THREADS + if ( sipSelfWasArg ) { + sipCpp->HYDROData_CalculationCase::SetStricklerTable( aRef ); + } else { + sipCpp->SetStricklerTable( aRef ); + } + Py_END_ALLOW_THREADS + } + %End /** * Add new one child region for calculation case. diff --git a/src/HYDROTools/CMakeLists.txt b/src/HYDROTools/CMakeLists.txt index f0515068..cceff83c 100644 --- a/src/HYDROTools/CMakeLists.txt +++ b/src/HYDROTools/CMakeLists.txt @@ -23,6 +23,7 @@ SET(PYFILES __init__.py controls.py interpolZ.py + interpolS.py ) # --- rules --- diff --git a/src/HYDROTools/interpolS.py b/src/HYDROTools/interpolS.py new file mode 100644 index 00000000..000df600 --- /dev/null +++ b/src/HYDROTools/interpolS.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- + +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 + +# Get current study id +salome.salome_init() +theStudy = salome.myStudy +theStudyId = salome.myStudyId + +"""Bad parameters exception""" +class BadParamsError(ValueError): + 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='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] + 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, medenum.MED_NO_DT, 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) \ No newline at end of file -- 2.39.2