docs: 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