def load_hours(network, min_load=0.9, max_load=1, boundaries=[0, 8760]):
    
    """Plot number of hours with line loading in selected range. 
    
    Parameters
    ----------
    network: PyPSA network container
        Holds topology of grid including results from powerflow analysis
    min_load: float
        Choose lower bound of relative load 
    max_load: float
        Choose upper bound of relative load
    boundaries: array
       Set boundaries of heatmap axis    
    
    """
    
    cmap_line = plt.cm.jet
    cmap_link = plt.cm.jet
    array_line = [['Line'] * len(network.lines), network.lines.index]

    load_lines = pd.Series(((abs(network.lines_t.p0[(
        abs(network.lines_t.p0.mul(network.snapshot_weightings, axis=0)) / 
        network.lines.s_nom_opt >= min_load) &
                                                    (
        abs(network.lines_t.p0.mul(network.snapshot_weightings, axis=0)) / 
        network.lines.s_nom_opt <= max_load)]) /
                            abs(network.lines_t.p0[(
                                abs(network.lines_t.p0) /
                                network.lines.s_nom_opt >= min_load) &
                                (abs(network.lines_t.p0) /
                                 network.lines.s_nom_opt <= max_load)]))
                            .sum()).data, index=array_line)

    array_link = [['Link'] * len(network.links), network.links.index]

    load_links = pd.Series(((abs(network.links_t.p0[(
        abs(network.links_t.p0.mul(network.snapshot_weightings, axis=0)) / 
        network.links.p_nom_opt >= min_load) &
                                                    (
        abs(network.links_t.p0.mul(network.snapshot_weightings, axis=0)) / 
        network.links.p_nom_opt <= max_load)]) /
                             abs(network.links_t.p0[(
                                abs(network.links_t.p0) /
                                network.links.p_nom_opt >= min_load) &
                                (abs(network.links_t.p0) /
                                 network.links.p_nom_opt <= max_load)]))
                            .sum()).data, index=array_link)

    load_hours = load_lines.append(load_links)

    ll = network.plot(
        line_colors=load_hours,
        line_cmap={
            'Line': cmap_line,
            'Link': cmap_link},
        bus_sizes=0,
        title="Number of hours with more then 90% load",
        line_widths=2)

    v1 = np.linspace(boundaries[0], boundaries[1], 101)
    v = np.linspace(boundaries[0], boundaries[1], 101)
    cb_Link = plt.colorbar(ll[2], boundaries=v1,
                           ticks=v[0:101:10])
    cb_Link.set_clim(vmin=boundaries[0], vmax=boundaries[1])

    cb = plt.colorbar(ll[1], boundaries=v,
                      ticks=v[0:101:10])
    cb.set_clim(vmin=boundaries[0], vmax=boundaries[1])

    cb.set_label('Number of hours')