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