docs: def storage_expansion(network, basemap=True, scaling=1, filename=None):
      """
      Plot storage distribution as circles on grid nodes
      Displays storage size and distribution in network.
      
      Parameters
      ----------
      network : PyPSA network container
          Holds topology of grid including results from powerflow analysis
      filename : str
          Specify filename
          If not given, figure will be show directly
      """
  
      stores = network.storage_units[network.storage_units.carrier == 
                                     'extendable_storage']
      batteries = stores[stores.max_hours == 6]
      hydrogen = stores[stores.max_hours == 168]
      storage_distribution =\
          network.storage_units.p_nom_opt[stores.index].groupby(
          network.storage_units.bus).sum().reindex(
                  network.buses.index, fill_value=0.)
      battery_distribution =\
          network.storage_units.p_nom_opt[batteries.index].groupby(
          network.storage_units.bus).sum().reindex(
                  network.buses.index, fill_value=0.)
      hydrogen_distribution =\
          network.storage_units.p_nom_opt[hydrogen.index].groupby(
          network.storage_units.bus).sum().reindex(
                  network.buses.index, fill_value=0.)
  
      sbatt = network.storage_units.index[
          (network.storage_units.p_nom_opt > 1) & (
                  network.storage_units.capital_cost > 10) & (
                      network.storage_units.max_hours == 6)]
      shydr = network.storage_units.index[
          (network.storage_units.p_nom_opt > 1) & (
                  network.storage_units.capital_cost > 10) & (
                      network.storage_units.max_hours == 168)]
  
      fig, ax = plt.subplots(1, 1)
      fig.set_size_inches(6, 6)
      
      msd_max = storage_distribution.max()
      msd_max_bat = battery_distribution.max()
      msd_max_hyd = hydrogen_distribution.max()
  
      if msd_max != 0:
          LabelVal = int(log10(msd_max))
      else:
          LabelVal = 0
      if LabelVal < 0:
          LabelUnit = 'kW'
          msd_max, msd_max_bat, msd_max_hyd = msd_max * \
              1000, msd_max_bat * 1000, msd_max_hyd * 1000
          battery_distribution = battery_distribution * 1000
          hydrogen_distribution = hydrogen_distribution * 1000
      elif LabelVal < 3:
          LabelUnit = 'MW'
      else:
          LabelUnit = 'GW'
          msd_max, msd_max_bat, msd_max_hyd = msd_max / \
              1000, msd_max_bat / 1000, msd_max_hyd / 1000
          battery_distribution = battery_distribution / 1000
          hydrogen_distribution = hydrogen_distribution / 1000
  
      if network.storage_units.p_nom_opt[sbatt].sum() < 1 and\
      network.storage_units.p_nom_opt[shydr].sum() < 1:
          print("No storage unit to plot")
      elif network.storage_units.p_nom_opt[sbatt].sum() > 1 and\
      network.storage_units.p_nom_opt[shydr].sum() < 1:
          network.plot(bus_sizes=battery_distribution * scaling, 
                       bus_colors='orangered', ax=ax, line_widths=0.3)
      elif network.storage_units.p_nom_opt[sbatt].sum() < 1 and\
      network.storage_units.p_nom_opt[shydr].sum() > 1:
          network.plot(bus_sizes=hydrogen_distribution * scaling,
                       bus_colors='teal', ax=ax, line_widths=0.3)
      else:
          network.plot(bus_sizes=battery_distribution * scaling,
                       bus_colors='orangered', ax=ax, line_widths=0.3)
          network.plot(bus_sizes=hydrogen_distribution * scaling,
                       bus_colors='teal', ax=ax, line_widths=0.3)
  
      if basemap and basemap_present:
          x = network.buses["x"]
          y = network.buses["y"]
          x1 = min(x)
          x2 = max(x)
          y1 = min(y)
          y2 = max(y)
  
          bmap = Basemap(resolution='l', epsg=network.srid, llcrnrlat=y1, urcrnrlat=y2, llcrnrlon=x1, urcrnrlon=x2, ax=ax)
          bmap.drawcountries()
          bmap.drawcoastlines()
  
      if msd_max_hyd !=0:
          plt.scatter([], [], c='teal', s=msd_max_hyd * scaling,
                  label='= ' + str(round(msd_max_hyd, 0)) + LabelUnit + ' hydrogen storage')
      if msd_max_bat !=0:
          plt.scatter([], [], c='orangered', s=msd_max_bat * scaling,
                  label='= ' + str(round(msd_max_bat, 0)) + LabelUnit + ' battery storage')
      plt.legend(scatterpoints=1, labelspacing=1, title='Storage size and technology', borderpad=1.3, loc=2)
      ax.set_title("Storage expansion")
  
      
      if filename is None:
          plt.show()
      else:
          plt.savefig(filename)
          plt.close()
  
      return