def stack_merge(s1, s2, out_stack_fn=None, sort=True, save=False):
    """Merge two stack objects
    """
    from pygeotools.lib import geolib
    from copy import deepcopy
    #Assumes input stacks have identical extent, resolution, and projection
    if s1.ma_stack.shape[1:3] != s2.ma_stack.shape[1:3]:
        print(s1.ma_stack.shape)
        print(s2.ma_stack.shape)
        sys.exit('Input stacks must have identical array dimensions')
    if not geolib.extent_compare(s1.extent, s2.extent):
        print(s1.extent)
        print(s2.extent)
        sys.exit('Input stacks must have identical extent')
    if not geolib.res_compare(s1.res, s2.res):
        print(s1.res)
        print(s2.res)
        sys.exit('Input stacks must have identical res')

    print("\nCombining fn_list and ma_stack")
    fn_list = np.array(s1.fn_list + s2.fn_list)

    if sort:
        #Sort based on filenames (should be datesort)
        sort_idx = np.argsort([os.path.split(x)[-1] for x in fn_list])
    else:
        sort_idx = Ellipsis

    #Now pull out final, sorted order
    fn_list = fn_list[sort_idx]
    ma_stack = np.ma.vstack((s1.ma_stack, s2.ma_stack))[sort_idx]
    #date_list = np.ma.dstack(s1.date_list, s2.date_list)
    #date_list_o = np.ma.dstack(s1.date_list_o, s2.date_list_o)
    source = np.array(s1.source + s2.source)[sort_idx]
    error = np.ma.concatenate([s1.error, s2.error])[sort_idx]
    #These are object arrays
    error_dict_list = np.concatenate([s1.error_dict_list, s2.error_dict_list])[sort_idx]

    print("Creating copy for new stack")
    s = deepcopy(s1)
    s.fn_list = list(fn_list)
    s.ma_stack = ma_stack
    s.source = list(source)
    s.error = error
    s.error_dict_list = error_dict_list
    #This will use original stack outdir
    if not out_stack_fn:
        s.get_stack_fn()
    else:
        s.stack_fn = out_stack_fn
    s.get_date_list()

    #These will preserve trend from one stack if present in only one stack
    #Useful when combining surface topo and bed topo
    if s1.datestack and s2.datestack:
        s.compute_dt_stats()
    if save and s1.datestack:
        s.write_datestack()

    if s1.stats and s2.stats:
        s.compute_stats()
    if save and s1.stats:
        s.write_stats()

    if s1.trend and s2.trend:
        s.compute_trend()
    if save and s1.trend:
        s.write_trend()

    if save:
        s.savestack()
    return s
Merge two stack objects