def get_sub_dim(src_ds, scale=None, maxdim=1024):
    """Compute dimensions of subsampled dataset 

    Parameters
    ----------
    ds : gdal.Dataset 
        Input GDAL Datset
    scale : int, optional
        Scaling factor
    maxdim : int, optional 
        Maximum dimension along either axis, in pixels
    
    Returns
    -------
    ns
        Numper of samples in subsampled output
    nl
        Numper of lines in subsampled output
    scale 
        Final scaling factor
    """
    ns = src_ds.RasterXSize
    nl = src_ds.RasterYSize
    maxdim = float(maxdim)
    if scale is None:
        scale_ns = ns / maxdim
        scale_nl = nl / maxdim
        scale = max(scale_ns, scale_nl)
    if scale > 1:
        ns = int(round(ns / scale))
        nl = int(round(nl / scale))
    return ns, nl, scale
