def calculate_mvgd_stats(nw):
    """
    MV Statistics for an arbitrary network

    Parameters
    ----------
    nw: :any:`list` of NetworkDing0
        The MV grid(s) to be studied

    Returns
    -------
    mvgd_stats : pandas.DataFrame
        Dataframe containing several statistical numbers about the MVGD
    """
    ##############################

    omega = 2 * pi * 50

    # close circuit breakers
    nw.control_circuit_breakers(mode='close')
    ##############################
    # Collect info from nw into dataframes
    # define dictionaries for collection
    trafos_dict = {}
    generators_dict = {}
    branches_dict = {}
    ring_dict = {}
    LA_dict = {}
    other_nodes_dict = {}
    lv_branches_dict = {}
    # initiate indexes
    trafos_idx = 0
    gen_idx = 0
    branch_idx = 0
    ring_idx = 0
    LA_idx = 0
    lv_branches_idx = 0

    # loop over MV grid districts
    for district in nw.mv_grid_districts():

        # node of MV station
        root = district.mv_grid.station()

        ###################################
        # get impedance of path to each terminal node
        # and get thermal capacity of first segment of path to each terminal node

        # store properties of terminal nodes in dictionaries
        # properties are e.g. impedance of path, length of path, thermal limit of first segment of path
        mv_impedances = {}
        mvlv_impedances = {}

        mv_path_lengths = {}
        mvlv_path_lengths = {}

        mv_thermal_limits = {}  # I_max of first segment on MV for each MV path
        lv_thermal_limits = {}  # I_max of first segment on LV for each LV path
        mvlv_thermal_limits = {}  # I_max of first segment on MV for each MVLV path

        n_outgoing_LV = 0
        n_stations_LV = 0

        n_outgoing_MV = 0

        G = district.mv_grid._graph

        for node in G.nodes():
            if isinstance(node, MVStationDing0):
                n_outgoing_MV += len(list(G.neighbors(node)))
                continue
            mv_impedance = 0
            mv_path_length = 0
            if not isinstance(node, MVCableDistributorDing0) and not isinstance(node, CircuitBreakerDing0):
                if not nx.has_path(G, root, node):
                    continue
                    #print(node, node.lv_load_area.is_aggregated) # only debug
                else:
                    path = nx.shortest_path(G, root, node)
                    for i in range(len(path) - 1):
                        mv_impedance += np.sqrt(
                            (G.adj[path[i]][path[i + 1]]['branch'].type[
                                 'L'] * 1e-3 * omega * \
                             G.adj[path[i]][path[i + 1]][
                                 'branch'].length) ** 2. + \
                            (G.adj[path[i]][path[i + 1]]['branch'].type[
                                 'R'] * \
                             G.adj[path[i]][path[i + 1]][
                                 'branch'].length) ** 2.)
                        mv_path_length += G.adj[path[i]][path[i + 1]][
                            'branch'].length

                    mv_impedances[node] = mv_impedance
                    mv_path_lengths[node] = mv_path_length
                    mv_thermal_limit = G.adj[path[0]][path[1]]['branch'].type['I_max_th']
                    mv_thermal_limits[node] = mv_thermal_limit

                    if isinstance(node, LVStationDing0):
                        # add impedance of transformers in LV station
                        lvstation_impedance = 0.
                        for trafo in node.transformers():
                            lvstation_impedance += 1. / np.hypot(trafo.r,trafo.x)  # transformers operating in parallel
                        if lvstation_impedance > 0.:  # avoid dividing by zero
                            lvstation_impedance = 1. / lvstation_impedance
                        else:
                            lvstation_impedance = 0.
                        # identify LV nodes belonging to LV station
                        for lv_LA in district.lv_load_areas():
                            for lv_dist in lv_LA.lv_grid_districts():
                                if lv_dist.lv_grid._station == node:
                                    G_lv = lv_dist.lv_grid._graph
                                    # loop over all LV terminal nodes belonging to LV station
                                    for lv_node in G_lv.nodes():
                                        if isinstance(lv_node, GeneratorDing0) or isinstance(lv_node, LVLoadDing0):
                                            path = nx.shortest_path(G_lv, node, lv_node)
                                            lv_impedance = lvstation_impedance
                                            lv_path_length = 0.
                                            for i in range(len(path)-1):
                                                lv_impedance += np.sqrt((G_lv.adj[path[i]][path[i+1]]['branch'].type['L'] * 1e-3 * omega * \
                                                                          G_lv.adj[path[i]][path[i+1]]['branch'].length)**2. + \
                                                                         (G_lv.adj[path[i]][path[i+1]]['branch'].type['R'] * \
                                                                          G_lv.adj[path[i]][path[i+1]]['branch'].length)**2.)
                                                lv_path_length += G_lv.adj[path[i]][path[i+1]]['branch'].length
                                            lv_thermal_limit = G_lv.adj[path[0]][path[1]]['branch'].type['I_max_th']

                                            mvlv_impedances[lv_node] = mv_impedance + lv_impedance
                                            mvlv_path_lengths[lv_node] = mv_path_length + lv_path_length
                                            lv_thermal_limits[lv_node] = lv_thermal_limit
                                            mvlv_thermal_limits[lv_node] = mv_thermal_limit

                                        elif isinstance(lv_node, LVStationDing0):
                                            n_outgoing_LV += len(list(G_lv.neighbors(lv_node)))
                                            n_stations_LV += 1

        # compute mean values by looping over terminal nodes
        sum_impedances = 0.
        sum_thermal_limits = 0.
        sum_path_lengths = 0.
        n_terminal_nodes_MV = 0

        # terminal nodes on MV
        for terminal_node in mv_impedances.keys():  # neglect LVStations here because already part of MVLV paths below
            if not isinstance(terminal_node, LVStationDing0) and not isinstance(terminal_node, MVStationDing0):
                sum_impedances += mv_impedances[terminal_node]
                sum_thermal_limits += mv_thermal_limits[terminal_node]
                sum_path_lengths += mv_path_lengths[terminal_node]
                n_terminal_nodes_MV += 1

        sum_thermal_limits_LV = 0.
        n_terminal_nodes_LV = 0

        # terminal nodes on LV
        for terminal_node in mvlv_impedances.keys():
            sum_impedances += mvlv_impedances[terminal_node]
            sum_thermal_limits += mvlv_thermal_limits[terminal_node]
            sum_thermal_limits_LV += lv_thermal_limits[terminal_node]
            sum_path_lengths += mvlv_path_lengths[terminal_node]
            n_terminal_nodes_LV += 1

        n_terminal_nodes = n_terminal_nodes_MV + n_terminal_nodes_LV

        if n_terminal_nodes < 1:
            mean_impedance = np.nan
            mean_thermal_limit = np.nan
            mean_path_length = np.nan
        else:
            mean_impedance = sum_impedances / n_terminal_nodes
            mean_thermal_limit = sum_thermal_limits / n_terminal_nodes
            mean_path_length = sum_path_lengths / n_terminal_nodes
        if n_terminal_nodes_LV < 1:
            mean_thermal_limit_LV = np.nan
        else:
            mean_thermal_limit_LV = sum_thermal_limits_LV / n_terminal_nodes_LV
        number_outgoing_LV = n_outgoing_LV  # / n_stations_LV
        number_outgoing_MV = n_outgoing_MV

        ###################################
        # compute path lengths (written by Miguel)
        max_mv_path = 0
        max_mvlv_path = 0

        # rings
        nodes_in_rings = []
        branches_in_rings = []
        for ring in district.mv_grid.rings_full_data():
            ring_idx += 1

            # generation cap
            ring_gen = 0
            for node in ring[2]:
                nodes_in_rings.append(node)
                if isinstance(node, GeneratorDing0):
                    ring_gen += node.capacity
            # length
            ring_length = 0
            for branch in ring[1]:
                branches_in_rings.append(branch)
                ring_length += branch.length / 1e3

            ring_dict[ring_idx] = {
                'grid_id': district.mv_grid.id_db,
                'ring_length': ring_length,
                'ring_capacity': ring_gen,
            }

        # transformers in main station
        for trafo in district.mv_grid.station().transformers():
            trafos_idx += 1
            trafos_dict[trafos_idx] = {
                'grid_id': district.mv_grid.id_db,
                's_max_a': trafo.s_max_a}

        # Generators and other MV special nodes
        cd_count = 0
        LVs_count = 0
        cb_count = 0
        lv_trafo_count = 0
        lv_trafo_cap = 0

        for node in district.mv_grid._graph.nodes():
            mv_path_length = 0
            mvlv_path_length = 0

            if isinstance(node, GeneratorDing0):
                gen_idx += 1
                isolation = not node in nodes_in_rings
                subtype = node.subtype
                if subtype == None:
                    subtype = 'other'
                generators_dict[gen_idx] = {
                    'grid_id': district.mv_grid.id_db,
                    'type': node.type,
                    'sub_type': node.type + '/' + subtype,
                    'gen_cap': node.capacity,
                    'v_level': node.v_level,
                    'isolation': isolation,
                }
                mv_path_length = district.mv_grid.graph_path_length(
                    node_source=root,
                    node_target=node)

            elif isinstance(node, MVCableDistributorDing0):
                cd_count += 1
            elif isinstance(node, LVStationDing0):
                LVs_count += 1
                lv_trafo_count += len([trafo for trafo in node.transformers()])
                lv_trafo_cap += np.sum([trafo.s_max_a for trafo in node.transformers()])

                if not node.lv_load_area.is_aggregated:
                    mv_path_length = district.mv_grid.graph_path_length(
                        node_source=root,
                        node_target=node)
                    max_lv_path = 0
                    for lv_LA in district.lv_load_areas():
                        for lv_dist in lv_LA.lv_grid_districts():
                            if lv_dist.lv_grid._station == node:
                                for lv_node in lv_dist.lv_grid._graph.nodes():
                                    lv_path_length = lv_dist.lv_grid.graph_path_length(
                                        node_source=node,
                                        node_target=lv_node)
                                    max_lv_path = max(max_lv_path, lv_path_length)
                    mvlv_path_length = mv_path_length + max_lv_path

            elif isinstance(node, CircuitBreakerDing0):
                cb_count += 1

            max_mv_path = max(max_mv_path, mv_path_length / 1000)
            max_mvlv_path = max(max_mvlv_path, mvlv_path_length / 1000)

        other_nodes_dict[district.mv_grid.id_db] = {
            'CD_count': cd_count,
            'LV_count': LVs_count,
            'CB_count': cb_count,
            'MVLV_trafo_count': lv_trafo_count,
            'MVLV_trafo_cap': lv_trafo_cap,
            'max_mv_path': max_mv_path,
            'max_mvlv_path': max_mvlv_path,
            'mean_impedance': mean_impedance,
            'mean_thermal_limit': mean_thermal_limit,
            'mean_thermal_limit_LV': mean_thermal_limit_LV,
            'mean_path_length': mean_path_length / 1.e3,
            'number_outgoing_LV': number_outgoing_LV,
            'number_outgoing_MV': number_outgoing_MV
        }

        # branches
        for branch in district.mv_grid.graph_edges():
            branch_idx += 1
            br_in_ring = branch['branch'] in branches_in_rings
            branches_dict[branch_idx] = {
                'grid_id': district.mv_grid.id_db,
                'length': branch['branch'].length / 1e3,
                'type_name': branch['branch'].type['name'],
                'type_kind': branch['branch'].kind,
                'in_ring': br_in_ring,
            }
        # Load Areas
        for LA in district.lv_load_areas():
            LA_idx += 1
            LA_dict[LA_idx] = {
                'grid_id': district.mv_grid.id_db,
                'is_agg': LA.is_aggregated,
                'is_sat': LA.is_satellite,
                # 'peak_gen':LA.peak_generation,
            }
            LA_pop = 0
            residential_peak_load = 0
            retail_peak_load = 0
            industrial_peak_load = 0
            agricultural_peak_load = 0
            lv_gen_level_6 = 0
            lv_gen_level_7 = 0
            for lv_district in LA.lv_grid_districts():
                LA_pop = + lv_district.population
                residential_peak_load += lv_district.peak_load_residential
                retail_peak_load += lv_district.peak_load_retail
                industrial_peak_load += lv_district.peak_load_industrial
                agricultural_peak_load += lv_district.peak_load_agricultural

                # generation capacity
                for g in lv_district.lv_grid.generators():
                    if g.v_level == 6:
                        lv_gen_level_6 += g.capacity
                    elif g.v_level == 7:
                        lv_gen_level_7 += g.capacity

                # branches lengths
                for br in lv_district.lv_grid.graph_edges():
                    lv_branches_idx += 1
                    lv_branches_dict[lv_branches_idx] = {
                        'grid_id': district.mv_grid.id_db,
                        'length': br['branch'].length / 1e3,
                        'type_name': br['branch'].type.to_frame().columns[0],  # why is it different as for MV grids?
                        'type_kind': br['branch'].kind,
                    }

            LA_dict[LA_idx].update({
                'population': LA_pop,
                'residential_peak_load': residential_peak_load,
                'retail_peak_load': retail_peak_load,
                'industrial_peak_load': industrial_peak_load,
                'agricultural_peak_load': agricultural_peak_load,
                'total_peak_load': residential_peak_load + retail_peak_load + \
                                   industrial_peak_load + agricultural_peak_load,
                'lv_generation': lv_gen_level_6 + lv_gen_level_7,
                'lv_gens_lvl_6': lv_gen_level_6,
                'lv_gens_lvl_7': lv_gen_level_7,
            })

        # geographic
        #  ETRS (equidistant) to WGS84 (conformal) projection
        proj = partial(
            pyproj.transform,
            # pyproj.Proj(init='epsg:3035'),  # source coordinate system
            # pyproj.Proj(init='epsg:4326'))  # destination coordinate system
            pyproj.Proj(init='epsg:4326'),  # source coordinate system
            pyproj.Proj(init='epsg:3035'))  # destination coordinate system
        district_geo = transform(proj, district.geo_data)
        other_nodes_dict[district.mv_grid.id_db].update({'Dist_area': district_geo.area})

    mvgd_stats = pd.DataFrame.from_dict({}, orient='index')
    ###################################
    # built dataframes from dictionaries
    trafos_df = pd.DataFrame.from_dict(trafos_dict, orient='index')
    generators_df = pd.DataFrame.from_dict(generators_dict, orient='index')
    other_nodes_df = pd.DataFrame.from_dict(other_nodes_dict, orient='index')
    branches_df = pd.DataFrame.from_dict(branches_dict, orient='index')
    lv_branches_df = pd.DataFrame.from_dict(lv_branches_dict, orient='index')
    ring_df = pd.DataFrame.from_dict(ring_dict, orient='index')
    LA_df = pd.DataFrame.from_dict(LA_dict, orient='index')

    ###################################
    # Aggregated data HV/MV Trafos
    if not trafos_df.empty:
        mvgd_stats = pd.concat([mvgd_stats, trafos_df.groupby('grid_id').count()['s_max_a']], axis=1)
        mvgd_stats = pd.concat([mvgd_stats, trafos_df.groupby('grid_id').sum()[['s_max_a']]], axis=1)
        mvgd_stats.columns = ['N° of HV/MV Trafos', 'Trafos HV/MV Acc s_max_a']

    ###################################
    # Aggregated data Generators
    if not generators_df.empty:
        # MV generation per sub_type
        mv_generation = generators_df.groupby(['grid_id', 'sub_type'])['gen_cap'].sum().to_frame().unstack(level=-1)
        mv_generation.columns = ['Gen. Cap. of MV ' + _[1] if isinstance(_, tuple) else _
                                 for _ in mv_generation.columns]
        mvgd_stats = pd.concat([mvgd_stats, mv_generation], axis=1)

        # MV generation at V levels
        mv_generation = generators_df.groupby(
            ['grid_id', 'v_level'])['gen_cap'].sum().to_frame().unstack(level=-1)
        mv_generation.columns = ['Gen. Cap. of MV at v_level ' + str(_[1])
                                 if isinstance(_, tuple) else _
                                 for _ in mv_generation.columns]
        mvgd_stats = pd.concat([mvgd_stats, mv_generation], axis=1)
        # Isolated generators
        mv_generation = generators_df[generators_df['isolation']].groupby(
            ['grid_id'])['gen_cap'].count().to_frame()  # .unstack(level=-1)
        mv_generation.columns = ['N° of isolated MV Generators']
        mvgd_stats = pd.concat([mvgd_stats, mv_generation], axis=1)

    ###################################
    # Aggregated data of other nodes
    if not other_nodes_df.empty:
        # print(other_nodes_df['CD_count'].to_frame())
        mvgd_stats['N° of Cable Distr'] = other_nodes_df['CD_count'].to_frame().astype(int)
        mvgd_stats['N° of LV Stations'] = other_nodes_df['LV_count'].to_frame().astype(int)
        mvgd_stats['N° of Circuit Breakers'] = other_nodes_df['CB_count'].to_frame().astype(int)
        mvgd_stats['District Area'] = other_nodes_df['Dist_area'].to_frame()
        mvgd_stats['N° of MV/LV Trafos'] = other_nodes_df['MVLV_trafo_count'].to_frame().astype(int)
        mvgd_stats['Trafos MV/LV Acc s_max_a'] = other_nodes_df['MVLV_trafo_cap'].to_frame()
        mvgd_stats['Length of MV max path'] = other_nodes_df['max_mv_path'].to_frame()
        mvgd_stats['Length of MVLV max path'] = other_nodes_df['max_mvlv_path'].to_frame()
        mvgd_stats['Impedance Z of path to terminal node (mean value)'] = \
            other_nodes_df['mean_impedance'].to_frame()
        mvgd_stats['I_max of first segment of path from MV station to terminal node (mean value)'] = \
            other_nodes_df['mean_thermal_limit'].to_frame()
        mvgd_stats['I_max of first segment of path from LV station to terminal node (mean value)'] = \
            other_nodes_df['mean_thermal_limit_LV'].to_frame()
        mvgd_stats['Length of path from MV station to terminal node (mean value)'] = \
            other_nodes_df['mean_path_length'].to_frame()
        mvgd_stats['Number of lines and cables going out from LV stations'] = \
            other_nodes_df['number_outgoing_LV'].to_frame()
        mvgd_stats['Number of lines and cables going out from MV stations'] = \
            other_nodes_df['number_outgoing_MV'].to_frame()

    ###################################
    # Aggregated data of MV Branches
    if not branches_df.empty:
        # km of underground cable
        branches_data = branches_df[branches_df['type_kind'] == 'cable'].groupby(
            ['grid_id'])['length'].sum().to_frame()
        branches_data.columns = ['Length of MV underground cables']
        mvgd_stats = pd.concat([mvgd_stats, branches_data], axis=1)

        # km of overhead lines
        branches_data = branches_df[branches_df['type_kind'] == 'line'].groupby(
            ['grid_id'])['length'].sum().to_frame()
        branches_data.columns = ['Length of MV overhead lines']
        mvgd_stats = pd.concat([mvgd_stats, branches_data], axis=1)

        # km of different wire types
        branches_data = branches_df.groupby(
            ['grid_id', 'type_name'])['length'].sum().to_frame().unstack(level=-1)
        branches_data.columns = ['Length of MV type ' + _[1] if isinstance(_, tuple) else _
                                 for _ in branches_data.columns]
        mvgd_stats = pd.concat([mvgd_stats, branches_data], axis=1)

        # branches not in ring
        total_br = branches_df.groupby(['grid_id'])['length'].count().to_frame()
        ring_br = branches_df[branches_df['in_ring']].groupby(
            ['grid_id'])['length'].count().to_frame()
        branches_data = total_br - ring_br
        total_br.columns = ['N° of MV branches']
        mvgd_stats = pd.concat([mvgd_stats, total_br], axis=1)
        branches_data.columns = ['N° of MV branches not in a ring']
        mvgd_stats = pd.concat([mvgd_stats, branches_data], axis=1)

    ###################################
    # Aggregated data of LV Branches
    if not lv_branches_df.empty:
        # km of underground cable
        lv_branches_data = lv_branches_df[lv_branches_df['type_kind'] == 'cable'].groupby(
            ['grid_id'])['length'].sum().to_frame()
        lv_branches_data.columns = ['Length of LV underground cables']
        mvgd_stats = pd.concat([mvgd_stats, lv_branches_data], axis=1)

        # km of overhead lines
        lv_branches_data = lv_branches_df[lv_branches_df['type_kind'] == 'line'].groupby(
            ['grid_id'])['length'].sum().to_frame()
        lv_branches_data.columns = ['Length of LV overhead lines']
        mvgd_stats = pd.concat([mvgd_stats, lv_branches_data], axis=1)

        # km of different wire types
        lv_branches_data = lv_branches_df.groupby(
            ['grid_id', 'type_name'])['length'].sum().to_frame().unstack(level=-1)
        lv_branches_data.columns = ['Length of LV type ' + _[1] if isinstance(_, tuple) else _
                                    for _ in lv_branches_data.columns]
        mvgd_stats = pd.concat([mvgd_stats, lv_branches_data], axis=1)

        # n° of branches
        total_lv_br = lv_branches_df.groupby(['grid_id'])['length'].count().to_frame()
        total_lv_br.columns = ['N° of LV branches']
        mvgd_stats = pd.concat([mvgd_stats, total_lv_br], axis=1)

    ###################################
    # Aggregated data of Rings
    if not ring_df.empty:
        # N° of rings
        ring_data = ring_df.groupby(['grid_id'])['grid_id'].count().to_frame()
        ring_data.columns = ['N° of MV Rings']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)

        # min,max,mean km of all rings
        ring_data = ring_df.groupby(['grid_id'])['ring_length'].min().to_frame()
        ring_data.columns = ['Length of MV Ring min']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)
        ring_data = ring_df.groupby(['grid_id'])['ring_length'].max().to_frame()
        ring_data.columns = ['Length of MV Ring max']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)
        ring_data = ring_df.groupby(['grid_id'])['ring_length'].mean().to_frame()
        ring_data.columns = ['Length of MV Ring mean']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)

        # km of all rings
        ring_data = ring_df.groupby(['grid_id'])['ring_length'].sum().to_frame()
        ring_data.columns = ['Length of MV Rings total']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)

        # km of non-ring
        non_ring_data = branches_df.groupby(['grid_id'])['length'].sum().to_frame()
        non_ring_data.columns = ['Length of MV Rings total']
        ring_data = non_ring_data - ring_data
        ring_data.columns = ['Length of MV Non-Rings total']
        mvgd_stats = pd.concat([mvgd_stats, ring_data.round(1).abs()], axis=1)

        # rings generation capacity
        ring_data = ring_df.groupby(['grid_id'])['ring_capacity'].sum().to_frame()
        ring_data.columns = ['Gen. Cap. Connected to MV Rings']
        mvgd_stats = pd.concat([mvgd_stats, ring_data], axis=1)
    ###################################
    # Aggregated data of Load Areas
    if not LA_df.empty:
        LA_data = LA_df.groupby(['grid_id'])['population'].count().to_frame()
        LA_data.columns = ['N° of Load Areas']

        mvgd_stats = pd.concat([mvgd_stats, LA_data], axis=1)

        LA_data = LA_df.groupby(['grid_id'])['population',
                                             'residential_peak_load',
                                             'retail_peak_load',
                                             'industrial_peak_load',
                                             'agricultural_peak_load',
                                             'total_peak_load',
                                             'lv_generation',
                                             'lv_gens_lvl_6',
                                             'lv_gens_lvl_7'
        ].sum()
        LA_data.columns = ['LA Total Population',
                           'LA Total LV Peak Load Residential',
                           'LA Total LV Peak Load Retail',
                           'LA Total LV Peak Load Industrial',
                           'LA Total LV Peak Load Agricultural',
                           'LA Total LV Peak Load total',
                           'LA Total LV Gen. Cap.',
                           'Gen. Cap. of LV at v_level 6',
                           'Gen. Cap. of LV at v_level 7',
                           ]
        mvgd_stats = pd.concat([mvgd_stats, LA_data], axis=1)

    ###################################
    # Aggregated data of Aggregated Load Areas
    if not LA_df.empty:
        agg_LA_data = LA_df[LA_df['is_agg']].groupby(
            ['grid_id'])['population'].count().to_frame()
        agg_LA_data.columns = ['N° of Load Areas - Aggregated']
        mvgd_stats = pd.concat([mvgd_stats, agg_LA_data], axis=1)

        sat_LA_data = LA_df[LA_df['is_sat']].groupby(
            ['grid_id'])['population'].count().to_frame()
        sat_LA_data.columns = ['N° of Load Areas - Satellite']
        mvgd_stats = pd.concat([mvgd_stats, sat_LA_data], axis=1)

        agg_LA_data = LA_df[LA_df['is_agg']].groupby(['grid_id'])['population',
                                                                  'lv_generation',
                                                                  'total_peak_load'].sum()
        agg_LA_data.columns = ['LA Aggregated Population',
                               'LA Aggregated LV Gen. Cap.', 'LA Aggregated LV Peak Load total'
                               ]
        mvgd_stats = pd.concat([mvgd_stats, agg_LA_data], axis=1)

    ###################################
    mvgd_stats = mvgd_stats.fillna(0)
    mvgd_stats = mvgd_stats[sorted(mvgd_stats.columns.tolist())]
    return mvgd_stats
MV Statistics for an arbitrary network

    Parameters
    ----------
    nw: :any:`list` of NetworkDing0
        The MV grid(s) to be studied

    Returns
    -------
    mvgd_stats : pandas.DataFrame
        Dataframe containing several statistical numbers about the MVGD