subcatchment.statistics not giving cumulative results

[anneliesesytsma]

Issue

I have been using the PySWMM subcatchment.statistics function for over a year now to get results from individual subcatchments, but recently it stopped providing cumulative (or reasonable) output. I'm running a very simply model consisting of (S1) impervious subcatchment routed to (S2) pervious subcatchment. I am interested in the % of runon to the pervious area that infiltrates. For this, I have been using S2.statistics to calculate the infiltration fraction, or infiltration/runon to the pervious area. However, when I run S2.statistics now, I get nearly 0 precipitation, infiltration, runon, etc. I get reasonable results when I use system_routing.runoff_stats, but this doesn't give me runon from impervious to the pervious subcatchment.

To Reproduce

with Simulation('inp_files/test.inp') as sim:
    S2 = Subcatchments(sim)["S2"]
    for step in sim:
        pass
    print(S2.statistics)

Expected behavior
With this code, I expect to get some values for evaporation, infiltration, runoff, precip, and runon. When I run this directly in SWMM, I get a total precip of 0.12 in, Total runon of 0.07 in, total infiltration of 0.19 in (see screen shot attached).

However, when I run the above code in PySWMM what I get is the following: {'evaporation': 0.0, 'infiltration': 0.0, 'peak_runoff_rate': 0.0, 'precipitation': 1.17877344e-315, 'runoff': 0.0, 'runon': 0.0}

I am attaching the .inp file (converted to .txt) I am using in case someone else is able to verify these results.

Desktop:
OS: Windows
Python Version: Python 3.6, 32-bit operating system
PySWMM Version 0.5.1

any help is appreciated!! Please let me know if you need more info/background.

test.txt

[bemcdonnell]

@anneliesesytsma, Sorry for the bug! Let me take a look at SWMM's code to see if maybe we mangled something in the latest release of SWMM. Two sets of eyes are better than one set if you'd like to investigate too

https://github.com/OpenWaterAnalytics/Stormwater-Management-Model/blob/854df626b32f74e5cfb3a07fbea3f51b5d7855bb/src/toolkitAPI.c#L2157-L2207

[katmratliff]

@bemcdonnell @anneliesesytsma I am looking at this too and have the same issue on my machine.

[katmratliff]

@anneliesesytsma do you recall the version of PySWMM (or the date downloaded/installed) you were last using when the subcatchment statistics was working as it should?

@bemcdonnell I wonder if something with how the SubcatchStats got refactored is causing problems with PySWMM if we did not update the wrapper accordingly, see here

[anneliesesytsma]

@katmratliff Unfortunately I don't recall..however I do update the PySWMM software pretty frequently (using !pip install) and the subcatchment statistics were working when I ran my code less than a month ago.

[barrc]

Maybe the addition of impervRunoff and pervRunoff to SubcatchStats in SWMM 5.1.13?

[anneliesesytsma]

@anneliesesytsma, Sorry for the bug! Let me take a look at SWMM's code to see if maybe we mangled something in the latest release of SWMM. Two sets of eyes are better than one set if you'd like to investigate too

https://github.com/OpenWaterAnalytics/Stormwater-Management-Model/blob/854df626b32f74e5cfb3a07fbea3f51b5d7855bb/src/toolkitAPI.c#L2157-L2207

@bemcdonnell I don't see anything that looks off in the code but I don't know what exactly to look for either..

[dickinsonre]

This is the way the stats for runoff are now defined in SWMM 5.1.013+ impervVol and PervVol have now been inserted in the list

// Input: j = subcatchment index
// rainVol = rainfall + snowfall volume (ft3)
// runonVol = runon volume from other subcatchments (ft3)
// evapVol = evaporation volume (ft3)
// infilVol = infiltration volume (ft3)
// impervVol = impervious runoff volume (ft3)
// pervVol = pervious runoff volume (ft3)
// runoffVol = runoff volume (ft3)
// runoff = runoff rate (cfs)

[bemcdonnell]

@barrc and @dickinsonre,

Thanks for the catch:

Yes SWMM 5.1.12 has the subcatchment statistics struct as:

//------------------------
// SUBCATCHMENT STATISTICS
//------------------------
typedef struct
{
    double       precip;
    double       runon;
    double       evap;
    double       infil;
    double       runoff;
    double       maxFlow;         
}  TSubcatchStats;

and 5.1.13 has it as:

//------------------------
// SUBCATCHMENT STATISTICS
//------------------------
typedef struct
{
    double       precip;
    double       runon;
    double       evap;
    double       infil;
    double       runoff;
    double       maxFlow;         
	double       impervRunoff;                                                 //(5.1.013)
	double       pervRunoff;                                                   //
}  TSubcatchStats;

Looks like we need to update SWMM and PySWMM.

SWMM should be updated in this function:

https://github.com/OpenWaterAnalytics/Stormwater-Management-Model/blob/854df626b32f74e5cfb3a07fbea3f51b5d7855bb/src/toolkitAPI.c#L2157-L2204

and PySWMM should be updated here:

https://github.com/OpenWaterAnalytics/pyswmm/blob/b35c50aeb855d665bb3c163892b1726dd54666ed/pyswmm/toolkitapi.py#L368-L379

and here:

https://github.com/OpenWaterAnalytics/pyswmm/blob/b35c50aeb855d665bb3c163892b1726dd54666ed/pyswmm/subcatchments.py#L626-L649

@anneliesesytsma, do you mind supporting these two projects to debug the issues? We will get to it eventually but if you need it sooner, we'd really appreciate your support!

[bemcdonnell]

@katmratliff, you also might be onto something with the double pointer. How are we handling this with the other stats functions?

[anneliesesytsma]

@bemcdonnell I'd be happy to but I admit that while I know how to edit the pyswmm code, I don't know how to edit the swmm code. where does the swmm code get edited?

[katmratliff]

@katmratliff, you also might be onto something with the double pointer. How are we handling this with the other stats functions?

The swmm_getSubcatchStats is the only stats function configured to use a double pointer. This is, however, how I've written the water quality info getters. So I think PySWMM does need to be updated accordingly to account for the double pointer.

[katmratliff]

@bemcdonnell it does seem more complicated than the water quality API though, and I'm not sure how to go about handling the double pointer with a structure.

[bemcdonnell]

@katmratliff, I wonder why we handled this stat different than the others..

https://github.com/OpenWaterAnalytics/Stormwater-Management-Model/blob/854df626b32f74e5cfb3a07fbea3f51b5d7855bb/src/stats.c#L793-L1023

Starting to think that we should comb through this and make some of these approaches more consistent

[katmratliff]

@bemcdonnell it wasn't originally written that way - it got changed around in this commit. My guess is that it got changed for SWIG but not sure...

[jennwuu]

@bemcdonnell it does seem more complicated than the water quality API though, and I'm not sure how to go about handling the double pointer with a structure.

@katmratliff

def subcatch_statistics(self, ID):
        """
        Get stats for a Subcatchment.

        :param str ID: Subcatchment ID
        :return: Group Stats
        :rtype: dict
        """
        index = self.getObjectIDIndex(tka.ObjectType.SUBCATCH.value, ID)

        # SWMM function handle.
        swmm_stats_func = self.SWMMlibobj.swmm_getSubcatchStats
        # SWMM function handle argument output structure.
        swmm_stats_func_arg = ctypes.POINTER(ctypes.POINTER(tka.SubcStats))
        # Define argument.
        swmm_stats_func.argtypes = (
            ctypes.c_int,
            swmm_stats_func_arg, )

        object_ptr = ctypes.POINTER(tka.SubcStats)()
        errcode = swmm_stats_func(
            ctypes.c_int(index), ctypes.byref(object_ptr))

        self._error_check(errcode)
        # Copy Items to Dictionary using Alias Names.
        out_dict = {}
        for attr in dir(object_ptr.contents):
            if "_" not in attr:
                out_dict[object_ptr.contents._py_alias_ids[attr]] = getattr(
                    object_ptr.contents, attr)

But I'm getting numbers from the rpt for runon and infiltration. I wonder if it is something weird with unit conversion?

int DLLEXPORT swmm_getSubcatchStats(int index, SM_SubcatchStats **subcatchStats)
///
/// Output:  Subcatchment Stats Structure (SM_SubcatchStats)
/// Return:  API Error
/// Purpose: Gets Subcatchment Stats and Converts Units
{
  int error_code_index = 0;
  double a;
  TSubcatchStats *tmp = (TSubcatchStats *)calloc(1, sizeof(TSubcatchStats));

  // Check if Open
  if (swmm_IsOpenFlag() == FALSE)
  {
    error_code_index = ERR_API_INPUTNOTOPEN;
  }

  // Check if Simulation is Running
  else if (swmm_IsStartedFlag() == FALSE)
  {
    error_code_index = ERR_API_SIM_NRUNNING;
  }

  // Check if object index is within bounds
  else if (index < 0 || index >= Nobjects[SUBCATCH])
  {
    error_code_index = ERR_API_OBJECT_INDEX;
  }
    // Copy Structure
  else
  {
    memcpy(tmp, stats_getSubcatchStat(index), sizeof(TSubcatchStats));
    *subcatchStats = (SM_SubcatchStats *)tmp;

    a = Subcatch[index].area;

    // Cumulative Runon Volume
    (*subcatchStats)->runon *= UCF(VOLUME);
    // Cumulative Infiltration Volume
    (*subcatchStats)->infil *= UCF(VOLUME);
    // Cumulative Runoff Volume
    (*subcatchStats)->runoff *= UCF(VOLUME);
    // Maximum Runoff Rate
    (*subcatchStats)->maxFlow *= UCF(FLOW);
    // Cumulative Rainfall Depth
    (*subcatchStats)->precip *= (UCF(RAINDEPTH) / a);
    // Cumulative Evaporation Volume
    (*subcatchStats)->evap *= UCF(VOLUME);
  }

    return error_getCode(error_code_index);
}

I think the volume calculations (such as evap, runon, runoff, etc.) is missing divide by subcatchment area and some unit conversion to get it into in/mm

void writeSubcatchRunoff()
{
    int    j;
    double a, x, r;

    if ( Nobjects[SUBCATCH] == 0 ) return;
    WRITE("");
    WRITE("***************************");
    WRITE("Subcatchment Runoff Summary");
    WRITE("***************************");
    WRITE("");
    fprintf(Frpt.file,

////////  Segment below modified for release 5.1.013.  /////////

"\n  ------------------------------------------------------------------------------------------------------------------------------"
"\n                            Total      Total      Total      Total     Imperv       Perv      Total       Total     Peak  Runoff"
"\n                           Precip      Runon       Evap      Infil     Runoff     Runoff     Runoff      Runoff   Runoff   Coeff");
    if ( UnitSystem == US ) fprintf(Frpt.file,
"\n  Subcatchment                 in         in         in         in         in         in         in    %8s      %3s",
        VolUnitsWords[UnitSystem], FlowUnitWords[FlowUnits]);
    else fprintf(Frpt.file,
"\n  Subcatchment                 mm         mm         mm         mm         mm         mm         mm    %8s      %3s",
        VolUnitsWords[UnitSystem], FlowUnitWords[FlowUnits]);
    fprintf(Frpt.file,
"\n  ------------------------------------------------------------------------------------------------------------------------------");

/////////////////////////////////////////////////////////////////

    for ( j = 0; j < Nobjects[SUBCATCH]; j++ )
    {
        a = Subcatch[j].area;
        if ( a == 0.0 ) continue;
        fprintf(Frpt.file, "\n  %-20s", Subcatch[j].ID);
        x = SubcatchStats[j].precip * UCF(RAINDEPTH);
        fprintf(Frpt.file, " %10.2f", x/a);
        x = SubcatchStats[j].runon * UCF(RAINDEPTH); 
        fprintf(Frpt.file, " %10.2f", x/a);
        x = SubcatchStats[j].evap * UCF(RAINDEPTH);
        fprintf(Frpt.file, " %10.2f", x/a);
        x = SubcatchStats[j].infil * UCF(RAINDEPTH); 
        fprintf(Frpt.file, " %10.2f", x/a);
        x = SubcatchStats[j].impervRunoff * UCF(RAINDEPTH);                    //(5.1.013)
        fprintf(Frpt.file, " %10.2f", x/a);                                    //
        x = SubcatchStats[j].pervRunoff * UCF(RAINDEPTH);                      //
        fprintf(Frpt.file, " %10.2f", x/a);                                    //
        x = SubcatchStats[j].runoff * UCF(RAINDEPTH);
        fprintf(Frpt.file, " %10.2f", x/a);
        x = SubcatchStats[j].runoff * Vcf;
        fprintf(Frpt.file, "%12.2f", x);
        x = SubcatchStats[j].maxFlow * UCF(FLOW);
        fprintf(Frpt.file, " %8.2f", x);
        r = SubcatchStats[j].precip + SubcatchStats[j].runon;
        if ( r > 0.0 ) r = SubcatchStats[j].runoff / r;
        fprintf(Frpt.file, "%8.3f", r);
    }
    WRITE("");
}

[bemcdonnell]

@jennwuu, values aside (whether volume or depth) I think it's fair to say that we're getting random memory garbage out of this function

[jennwuu]

@bemcdonnell I think if you convert the python code to use double pointer for the struct (to match the api) and fix the unit conversion in the SWMM toolkitapi function, then it should work

[bemcdonnell]

@jennwuu, if you fix it I will help pay for food for your sloths :-)

[anneliesesytsma]

@jennwuu @bemcdonnell I will also help pay for sloth food if you fix this!!

[jennwuu]

I will leave the unit conversion alone so we don't have to update the swmm.dll @bemcdonnell @katmratliff @anneliesesytsma

This means runoff, runon, and infiltration will be in volume [ft3/m3] instead of volume depth [in/mm].

Check out PR #212.

[katmratliff]

@jennwuu, values aside (whether volume or depth) I think it's fair to say that we're getting random memory garbage out of this function

@bemcdonnell @jennwuu I was getting similar garbage values when I was testing and trying to figure out how to deal with the double pointers in the rest of the WQ API.

[bemcdonnell]

@katmratliff, did it get resolved?

[jennwuu]

@jennwuu, values aside (whether volume or depth) I think it's fair to say that we're getting random memory garbage out of this function

@bemcdonnell @jennwuu I was getting similar garbage values when I was testing and trying to figure out how to deal with the double pointers in the rest of the WQ API.

@katmratliff

I think you deal with double pointers by passing the address of ctype.pointer(structure).
https://github.com/OpenWaterAnalytics/pyswmm/blob/24651d377832198b895ee0fa41ecfe23c0d0c954/pyswmm/swmm5.py#L1805-L1807

The arg of double pointer should ctype.pointer(ctype.pointer(structure))
https://github.com/OpenWaterAnalytics/pyswmm/blob/24651d377832198b895ee0fa41ecfe23c0d0c954/pyswmm/swmm5.py#L1799

[katmratliff]

@bemcdonnell @jennwuu yes and yes! It's not a problem now, was just recalling that the values I was getting while trying to get the double pointers working last year (before things were working) were similar to the values that @anneliesesytsma was getting before the fix (zeros and tiny numbers). All working now, sorry for the confusion there!

[bemcdonnell]

@anneliesesytsma, Looks like @jennwuu has just released v0.5.2 to PyPI. Let us know if you have any other problems

[anneliesesytsma]

@bemcdonnell @katmratliff @jennwuu THANK YOU!! It works!!!