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