HEC-RAS批处理的实现

简要说明

利用MATLAB等软件控制HEC-RAS的输入输出操作可以实现批处理功能，大幅节约人为操作时间。今晚整理了几篇批处理的参考文献，如下所示。

• 文献【1】的2.3.3小节有一个详细的方法介绍；
• 文献【2，3】里有MATLAB代码有下载网址，但现在貌似要60美刀。
• 文献【4】给出了python脚本。
• 文献【5】给出了结合BP算法（python）的二维模型糙率自动率定方法。

1.借助pywin32库查询得到HEC-RAS的COM模块；2.将相应的COM模块导入Python环境,对HEC-RAS的调用、改写和运算；3.通过Python中h5py库,对HEC-RAS二维模型计算结果进行调用和读取。

因此可以参照文献【1】和【4】的思路完成代码设计，实现HEC_RAS的MATLAB操作。

注意：利用MATLAB的actxserver命令创建 COM 服务器。记住这里要确认软件COM组件安装成功，并开放了许可。

MATLAB代码

%==============================================================================
% FluEgg -Fluvial Egg Drift Simulator
%==============================================================================
% Copyright (c) 2013 Tatiana Garcia

   % This program is free software: you can redistribute it and/or modify
    % it under the terms of the GNU General Public License version 3 as published by
    % the Free Software Foundation (currently at http://www.gnu.org/licenses/agpl.html) 
    % with a permitted obligation to maintain Appropriate Legal Notices.

    % This program is distributed in the hope that it will be useful,
    % but WITHOUT ANY WARRANTY; without even the implied warranty of
    % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    % GNU General Public License for more details.

    % You should have received a copy of the GNU General Public License
    % along with this program.  If not, see <http://www.gnu.org/licenses/>.

%%%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::%
%%                       Extract output from HEC-RAS                      %
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::%
%-------------------------------------------------------------------------%
% This function is used to extract data from a HEC-RAS project and produces %
% References:                                                             %
% Goodell, C.R. 2014.                                                     %
% Breaking the HEC-RAS Code: A User's Guide to Automating HEC-RAS. A User's
% Guide to Automating HEC-RAS. h2ls. Portland, OR.                        %
%-------------------------------------------------------------------------%
%                                                                         %
%-------------------------------------------------------------------------%
%   Created by      : Tatiana Garcia                                      %
%   Date            : March 29, 2016                                      %
%   Last Modified   : April 18, 2016                                      %
%-------------------------------------------------------------------------%
% Inputs:
% Outputs:
%  
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::%

function [Riverinputfile]=Extract_RAS(strRASProject,handles,profile)

%% Creates a COM Server for the HEC-RAS Controller

% The command above depends on the version of HEC-RAS you have, in my case
% I am using version 5.0.3
try
    RC = actxserver('RAS503.HECRASController');
catch
    try %HECRAS 5.0.2
        RC = actxserver('RAS502.HECRASController');
    catch
        try %HECRAS 5.0.1
            RC = actxserver('RAS501.HECRASController');
        catch
            try %HECRAS 5.0.0
                RC = actxserver('RAS500.HECRASController');
            catch
            end %HECRAS 5.0.0
        end %HECRAS 5.0.1
    end%HECRAS 5.0.2
end %HECRAS 5.0.3


%% Open the project
%strRASProject = 'D:\Asian Carp\Asian Carp_USGS_Project\Tributaries data\Sandusky River\SANDUSKY_Hec_RAS_mod\Sandusky_mod_II\BallvilleDam_Updated.prj';
RC.Project_Open(strRASProject); %open and show interface, no need to use RC.ShowRAS in Matlab

%% Define Variables
% get variables from GUI with user selection
lngRiverID = get(handles.popup_River,'Value');   % RiverID
lngReachID = get(handles.popup_Reach,'Value');   % ReachID
%%
if profile==-1
lngProfile = get(handles.popup_HECRAS_profile,'Value')-1;   % Profile Number, 1 is all profiles
else
    %The first profile in HEC-RAS corresponds to MAX WS, then the profile
    % corresponding to initial conditions is:
    lngProfile=profile+1;    
end    

lngUpDn = 0;      % Up/Down index for nodes with multiple sections (only used for bridges)
% Output ID of Variables ( see page 247 in reference book for more details)
% lngWS_ID = 2;                     % The Water Surface Elevation ID is 2.
lngVelChnl_ID = 25;                 % The  ID of the average velocity of flow for the main channel is 23.
lngHydrDepthC_ID = 128;             % Hydraulic depth in channel.
lngQChannel_ID = 7;                 % Flow in the main channel.
lngHydrRadiusC_ID = 210;            % ID for hydraulic radius in channel.
lngMannWtdChnl_ID = 45;             % ID for Conveyance weighted Manning's n for the main channel.
lngLengthChnl_ID = 42;              % ID for Downstream reach length of the main channel to next XS
                                    %(unless BR in d/s, then this is the distance to the deck/roadway).

lngNum_XS = RC.Schematic_XSCount(); %Number of XS - HEC-RAS Controller will populate.

[~,~,lngNum_RS]=RC.Geometry_GetNodes(lngRiverID,lngReachID,0,0,0);%Number of nodes
        
% Preallocate memory
strRS = cell(lngNum_RS,1); %Array of names of the nodes-->River station name.  See page 36 in book
strNodeType = strRS;       %Pre-allocate array for node type

% Preallocate memory for output vectors
%sngWS = nan(lngNum_RS,1);
sngVelChnl = nan(lngNum_XS,1);
sngHydrDepthC = nan(lngNum_XS,1);
sngQChannel = nan(lngNum_XS,1);
sngHydrRadiusC = nan(lngNum_XS,1);
sngMannWtdChnl = nan(lngNum_XS,1);
sngLengthChnl = nan(lngNum_XS,1);

%% Run the current plan
%RC.Compute_CurrentPlan(0,0); %from book: = RC.Compute_CurrentPlan(lngMessages,strMessages(), True)
% RC.Compute_HideComputationWindow; %To hide Computation Window

% Uncoment the lines below for info about current plan and project
% Do not delete, we might need this in the future
%==========================================================================
%Current_Plan= RC.CurrentPlanFile()
%Current_Geometry_File= RC.CurrentGeomFile()
%[lngNumProf,strProfileName]=RC.Output_GetProfiles(0,0) %Profiles names and todatl number

%% Output Results
XC_counter=0;
% Extracts variable info from XS to Xs
for i=1:lngNum_RS
    strNodeType{i}=RC.Geometry.NodeCType(lngRiverID,lngReachID,i);
    
    if strcmp(strNodeType{i},'')% An empty strings (i.e. ' ') denotes a cross section;
            XC_counter=XC_counter+1;
        % Here we assing the river station name to each node
        strRS{XC_counter} = RC.Geometry.NodeRS(lngRiverID,lngReachID,i);
        % Extracts average velocity of flow for the main channel
        sngVelChnl(XC_counter) = abs(RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngVelChnl_ID));
        % Extracts hydraulic depth in channel
        sngHydrDepthC(XC_counter) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngHydrDepthC_ID);
        % Extracts flow in the main channel.
        sngQChannel(XC_counter) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngQChannel_ID);
        % Extracts hydraulic radius in channel.
        sngHydrRadiusC(XC_counter) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngHydrRadiusC_ID);
        % Extracts conveyance weighted Manning's n for the main channel.
        sngMannWtdChnl(XC_counter) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngMannWtdChnl_ID);
        % Extracts conveyance weighted Manning's n for the main channel.
        sngLengthChnl(XC_counter) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
            lngProfile,lngLengthChnl_ID);       
        % Extracts Water Surface Elevation at each XS
        %sngWS(i) = RC.Output_NodeOutput(lngRiverID,lngReachID,i,lngUpDn,...
        %           lngProfile,lngWS_ID);       
    end%if is a cross section
    
end %for
sngLengthChnl(end)=0; %There is not data of length channel for the last cell

%% Generate Rivirinputfile dataset for the FluEgg model
Riverinputfile=Generate_Rivirinputfile();

try
    %Kill Hec-ras 
    !taskkill /im ras.exe
    RC.QuitRAS
    
catch
end
delete(RC);



%% :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    function Riverinputfile=Generate_Rivirinputfile()
        CumlDistance=[0;cumsum(sngLengthChnl(1:end-1))/1000];%In km
        CumlDistance=[((CumlDistance(1:end-1)+CumlDistance(2:end))/2); CumlDistance(end)]; %FluEgg cells are located between 2 XS
        CellNumber=(1:1:lngNum_XS)';
        ks=(8.1.*sngMannWtdChnl.*sqrt(9.81)).^6;
        Ustar=abs(sngVelChnl./(8.1*((sngHydrRadiusC./ks).^(1/6))));
        %% Temperature choice
        Temperature=str2double(get(handles.Const_Temp(2),'String'));
        if isnan(Temperature) %If user didn't input temperature
            ed = errordlg('Please input temperature','Error');
            set(ed, 'WindowStyle', 'modal');
            uiwait(ed);
            return
        end
        Riverinputfile=[CellNumber CumlDistance sngHydrDepthC sngQChannel...
                        sngVelChnl zeros(lngNum_XS,1) zeros(lngNum_XS,1) Ustar ...
                        Temperature*ones(lngNum_XS,1)]; % default temperature 22C
    end
%% :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

end
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::%
%% <<<<<<<<<<<<<<<<<<<<<<<<< END OF FUNCTION >>>>>>>>>>>>>>>>>>>>>>>>>>>>%%
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::%

%% Notes:
% When there is not data the controller puts a very large number: 3.39999995214436e+38

挖个坑，急用可优先考虑python，现在寒潮。。。冻。

了了于心，如如不动。

参考文献

[1]. 冯诗韵,王飞儿,俞洁.基于Matlab软件自动化求取参数的HEC-RAS模型构建[J].环境科学学报,2020,40(02):623-630.
[2]. Leon A S , Goodell C . Controlling HEC-RAS using MATLAB[J]. Environmental Modelling & Software, 2016, 84(OCT.):339-348.
[3]. Goodell C . Breaking the HEC-RAS Code - A User’s Guide to Automating HEC-RAS[M]. 2014.
[4]. Dysarz T . Application of Python Scripting Techniques for Control and Automation of HEC-RAS Simulations[J]. Water, 2018, 10(10):1382.
[5].夏铭辉,秦景,牛文龙,雷添杰.基于BP神经网络的HEC-RAS二维模型糙率参数自动率定[J].水利水电技术,2020,51(05):38-46.