Extract statistics from .res11 or .res1d-files

Here is a code example on how to extract statistics from a res11 or res1d-files. Copy the code and the replace the file names and paths.

# -*- coding: utf-8 -*-
import sys
import os
from hydroinform.ts_tools import HydroTs
from hydroinform import DFS, Mike_river
import pandas as pd
import numpy as np


def get_single_point_ts(filename, rivername, topo_id, chainage):
    """
    Gets a single timeseries from a single point in a river
    Note, this method reads the entire output file, so use with care
    """
    with Mike_river.River_setup(filename) as res:
        reach = next((r for r in res.reaches if r.name == 'rivername' and r.topo_id==topo_id), None)
        if reach != None:
            gp = next((p for p in reach.grid_points if p.chainage == chainage), None)
            ts = res.get_timeseries_Q(gp)
            return ts

def statistics(filename, prename):
    """
    Returns a dataframe with statistics on the water flow from either a .res11- or res1d-file
    """
    res11=Mike_river.River_setup.from_file(filename)

    columns=['Name', 'TopoId', 'Branch', 'Downstream','Chainage','X','Y', 'Min', 'Max', 'Average', 'DFFV_EQR', 'DVFI_EQR', 'DVPI_EQR', 'Q95', 'Q75', 'Q50', 'Q25', 'MedianMin']    
    nbqpoints = np.sum([len(r.grid_points_q) for r in res11.reaches])
    data = pd.DataFrame(columns=columns,index = range(nbqpoints))
   
    index =0
    for reach in res11.reaches: #Loop all reaches
        for gp in reach.grid_points_q: #Loop all q-points
            ts = res11.get_timeseries_Q(gp) # get the discharge time series
            hts = HydroTs(ts) # Create a hydro time series
            #set all values
            idstring=prename + reach.topo_id + '_' + reach.name +'_' + str(round(gp.chainage,1))
            if(gp.downstream!=''):
                idstring = idstring +'_' + gp.downstream

            data.at[index, 'Name'] = prename + gp.name
            data.at[index, 'TopoId'] = reach.topo_id
            data.at[index, 'Branch']=reach.name
            data.at[index, 'Chainage']=gp.chainage
            data.at[index, 'Downstream']=gp.downstream
            data.at[index, 'X']=gp.x
            data.at[index, 'Y']=gp.y
            data.at[index, 'Min']=ts.min()
            data.at[index, 'Max']=ts.max()
            data.at[index, 'Average']=np.average(ts)
            data.at[index, 'DFFV_EQR']=hts.DFFV_EQR()
            data.at[index, 'DVFI_EQR']=hts.DVFIEQR()
            data.at[index, 'DVPI_EQR']=hts.DVPIEQR()
            data.at[index, 'Q95']=hts.Q95
            data.at[index, 'Q75']=hts.Q75
            data.at[index, 'Q50']=hts.Q50
            data.at[index, 'Q25']=hts.Q25
            data.at[index, 'MedianMin']=hts.medianmin()
            index+=1
    res11.dispose()
    return data



if __name__ == "__main__":

    #get_single_point_ts(r'QPoints2\R201803_vandweb\m11_res\abs\DK1_Abs.res11','MERNAA', 'novana_model', 8251).to_csv(r'c:\temp\testData.csv')

    import shapefile
    df =pd.DataFrame()

    MikeHydroFiles=[
           r'C:\projects\GEUS\DKM2019\production\DK1\DK1_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK2\DK2_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK3\DK3_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK4\DK4_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK5\DK5_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK6\DK6_mh_2018.res1d',
           r'C:\projects\GEUS\DKM2019\production\DK7\DK7_mh_2018.res1d'
           ]

    i=0
    #Loop the files
    for file in MikeHydroFiles:
        print('Reading: ' + file)
        i=i+1
        temp=statistics(file, 'DKM'+str(i)+'_')
        temp['DKModel']=i
        #Join the dataframes
        df=df.append(temp, ignore_index=True)
    
    #Write to shapefile
    with shapefile.Writer(r'c:\temp\output.shp', shapeType=shapefile.POINT) as w:
        w.autoBalance=True
        #First columns are text (Branch name)
        for c in df.columns[0:4]:
            w.field(c, 'C')
        #Remaining columns are numbers
        for c in df.columns[4:-1]:
            w.field(c, 'F', decimal=10)
            
        #Dkmodel number is integer    
        w.field(df.columns[-1], 'N')

        for i, row in df.iterrows():
            w.record(*row)
            w.point(row['X'],row['Y'])