Script tool: Create storage area footprints

This script tool takes the maximum surface area of each storage unit, generates a circular polygon of equivalent surface area, and creates a background layer in the Map panel depicting the approximate maximum aerial extent, or 'footprint,' of each storage. This tool is designed to work with PCSWMM version 7.6 or later.

An example output from the Create storage area footprints tool:

How the script works

This script uses the buffer tool to generate a circular polygon that is of approximately the same area as the maximum surface area specified by each storage curve.

For storages with TABULAR storage curves, this maximum value is the maximum area specified in the Curves Editor.

For storages with FUNCTIONAL storage curves, this maximum value is the area value corresponding to the full depth of the storage, defined by the depth-area relationship of the storage.

The generated footprint layer contains the following attributes:

Limitations and caveats

Because the buffer tool relies on the units of the Map Coordinate System to calculate distances, the tool will show an error if the coordinate system is in units of degree, or if the coordinate system units are inconsistent with the simulation units. For example, if the simulation units are in CFS (Cubic Feet per Second) but the coordinate system is in meters, the script will prompt the user to reproject the layer to a new coordinate system before running the tool.

Discrepancies between the MAXAREA attribute of the footprint and calculated GIS area attributes may be due to projection approximations defined by the selected Coordinate System, or due to limitations of the polygon buffer, which approximates a circle by generating a polygon with many edges.

How to use this script

This tool is designed to work with PCSWMM version 7.6 or later.

The user interface of this script is shown below.

1. Select the output file location and name of the footprint layer. By default, the tool will save the footprint layer to a shapefile called "[Model Name].StorageFootprints" in the model folder when 'Create' is clicked, unless the user changes the output file name in the textbox.
2. Click 'Selected entities only' to create footprints for only the selected storages. Clicking 'Create' will create a footprints shapefile in the output location, and a layer in the Map panel of the project.

Download script

Click the Copy button to copy the script and add it to the Script Editor in PCSWMM.

# tags: MapToolsSWMM,Article
# created on: 2022-04-28
# updated on: 2023-05-24
# name: Create storage area footprints
# file: create_storage_area_footprints
# description: Create circular footprints representing the maximum surface area of storages

import os
import math

# ----------------------------------------------------------
# global functions
# ----------------------------------------------------------
def add_layer(out_lyr, shape):
    """ create and new layer in the map """ 
    # if the layer already exists, return it
    for lyr in pcpy.Map.Layers:
        if lyr.FilePath == out_lyr:
            return lyr
    # if the file exists in the disk, open it
    if os.path.isfile(out_lyr):
        return pcpy.Map.open_layer(out_lyr)
    # if the file does not exist, create a new one
    return pcpy.Map.add_layer(out_lyr, shape)

def initialize_footprint_layer(layer):
    """ generate and format and new layer for storage footprints """
    layer.delete_entities() #clear so that previous iterations of the tool run are deleted/overwritten   
    attributes = layer.Attributes
    if 'NAME' not in attributes:
        layer.add_attribute(pcpy.Attribute('NAME','Text'))
    if 'DEPTH' not in attributes:
        layer.add_attribute(pcpy.Attribute('DEPTH','Number'))
    if 'MAXDEPTH' not in attributes:
        layer.add_attribute(pcpy.Attribute('MAXDEPTH','Number'))
    if 'MAXAREA' not in attributes:
        layer.add_attribute(pcpy.Attribute('MAXAREA','Number'))

    layer.Attributes['DEPTH'].Units = '<ModelDepth>'
    layer.Attributes['MAXDEPTH'].Units = '<ModelDepth>'
    layer.Attributes['MAXAREA'].Units = '<ModelHead>\xB2'


# ----------------------------------------------------------
# Footprint 
# ----------------------------------------------------------
class Footprint():
    """
    create a circular footprint representing the maximum area of the storage curve, 
    for both TABULAR and FUNCTIONAL storage curves, and add it to the generated 'Storage Footprint' layer.
    """   
    def __init__(self, src_entity, curves):
        self.sto_entity = src_entity
        self.curves = curves
      
    def get_max_curve_area(self):
        """ get maximum area value from the storage curve""" 
        if self.sto_entity['ShapeCurve'] == 'FUNCTIONAL':
            A = self.sto_entity['Coeff']
            B = self.sto_entity['Exponent']
            C = self.sto_entity['Constant']
            return (A * (self.sto_entity['Depth']**B) + C)
            
        else:
            curve_name = self.sto_entity['CurveName']
            curve = self.curves.get(curve_name)
            if curve is None:
                return 0       
            data = curve.Data
            areas = [row.Y for row in data]
            return max(areas)

    def generate_storage_footprint(self, max_curve_area):
        """ buffer the geometry of the storage to represent the maximum curve area """ 
        radius = math.sqrt(max_curve_area/math.pi)
        # pull the geometry from entities layer for matching storage name:
        sto_geom = self.sto_entity.Geometry
        centr = sto_geom.Centroid # to handle both point and polygon storage geometries
        pygeom_center = pcpy.Point(centr.X, centr.Y)
     
        buffered_geom = pygeom_center.buffer(radius)
        return buffered_geom

    def to_layer(self, layer):
        """ add footprint entities to a layer and format the entities """
        max_curve_area = self.get_max_curve_area()
        if max_curve_area == 0:
            print('Cannot calculate max curve area for storage "%s"' % self.sto_entity['Name'])
            return
        buffered_geom = self.generate_storage_footprint(max_curve_area)    
        foot_entity = layer.add_entity(buffered_geom)
        foot_entity['NAME'] = self.sto_entity['Name']
        foot_entity['DEPTH'] = self.sto_entity['Depth']
        foot_entity['MAXDEPTH'] = self.sto_entity['MaxDepth']
        foot_entity['MAXAREA'] = max_curve_area
        
# ----------------------------------------------------------
# Callbacks
# ----------------------------------------------------------
def on_button_click(sender, e):
    curves = pcpy.SWMM.Curves       
    try:
        out_lyr = shp_fname.Text      
        if not out_lyr.endswith('.shp'):
            out_lyr = out_lyr + '.shp'
            
        new_layer = add_layer(out_lyr,'Polygon')
        initialize_footprint_layer(new_layer)
        
        if ((flow_units in ['CFS','GPM','MGD']) and (storage_layer.CoordinateSystem.Units != 'ft')):
            raise IOError("This model uses simulation units that are inconsistent with the units of the Storages layer coordinate system (%s).\n" \
                     "Reproject the layer to a Projected CRS with units of feet (ft) for accurate footprints." %storage_layer.CoordinateSystem.Units)   
        if ((flow_units in ['CMS','LPS','MLD']) and (storage_layer.CoordinateSystem.Units != 'm')):
            raise IOError("This model uses simulation units that are inconsistent with the units of the Storages layer coordinate system (%s).\n" \
                     "Reproject the layer to a Projected CRS with units of meters (m) for accurate footprints." %storage_layer.CoordinateSystem.Units)   

        src_entities = storage_layer.get_selected_entities() if chk_selected_only.Checked else storage_layer.get_entities()
        for src_entity in src_entities:
            footprint = Footprint(src_entity, curves)
            footprint.to_layer(new_layer)
               
        pcpy.Map.refresh()       
        fm.close()
    except Exception as err:
        pcpy.show_messagebox(str(err), '', pcpy.Enum.IconType.Error)

def on_form_load(sender, e):
    # set checked only status:
    has_selection = len(storage_layer.get_selected_entities()) > 0
    chk_selected_only.Enabled = has_selection
    chk_selected_only.Checked = has_selection
    
# ---------------- Main ---------------------
pcpy.clear_output()
storage_layer = pcpy.Map.Layer['Storages']
flow_units = pcpy.SWMM.Options.FlowUnits

fm = pcpy.Form('Create Storage Area Footprints')
fm.HelpLink='https://support.chiwater.com/171657'
fm.add_title('Create Storage Area Footprints',
"""
Create a circular polygon 'footprint' depicting the maximum surface area of each storage unit.
""")
fm.BoxWidth = 500
out_path = os.path.join('%s.StorageFootprints.shp' % pcpy.SWMM.FilePath[0:-4])
shp_fname = fm.add_textbox_with_button('Output file','save',out_path)
chk_selected_only = fm.add_checkbox('Selected storage nodes only')
fm.Button1.Text = 'Create'
fm.Button1.Click += on_button_click
fm.Form.Load += on_form_load
fm.show()

Disclaimer

All scripts are provided as-is, and no warranty is made as to the script's accuracy, completeness, or suitability for any particular purpose. This script may contain errors and require correction or modification by the modeler before its use.

The modeler is obligated to review the script code in its entirety and determine its suitability for its intended usage.

This script may also be updated from time to time to fix issues or provide enhanced functionality. If you have suggestions for improvement, please create a support ticket to let us know.