a function to setup using all values from a
    layer-stress period pair for observations.  Writes
    an instruction file and a _setup_ csv used
    construct a control file.

    Parameters
    ----------
    hds_file : str
        a MODFLOW head-save file.  If the hds_file endswith 'ucn',
        then the file is treated as a UcnFile type.
    kperk_pairs : iterable
        an iterable of pairs of kper (zero-based stress
        period index) and k (zero-based layer index) to
        setup observations for.  If None, then a shit-ton
        of observations may be produced!
    skip : variable
        a value or function used to determine which values
        to skip when setting up observations.  If np.scalar(skip)
        is True, then values equal to skip will not be used.
        If skip can also be a np.ndarry with dimensions equal to the model.
        Obscervations are set up only for cells with Non-zero values in the array.
        If not np.ndarray or np.scalar(skip), then skip will be treated as a lambda function that
        returns np.NaN if the value should be skipped.
    prefix : str
        the prefix to use for the observation names. default is "hds".

    Returns
    -------
    (forward_run_line, df) : str, pd.DataFrame
        a python code str to add to the forward run script and the setup info for the observations

    Note
    ----
    requires flopy

    writes <hds_file>.dat.ins instruction file

    writes _setup_<hds_file>.csv which contains much
    useful information for construction a control file

def setup_hds_obs(hds_file,kperk_pairs=None,skip=None,prefix="hds"):
    """a function to setup using all values from a
    layer-stress period pair for observations.  Writes
    an instruction file and a _setup_ csv used
    construct a control file.

    Parameters
    ----------
    hds_file : str
        a MODFLOW head-save file.  If the hds_file endswith 'ucn',
        then the file is treated as a UcnFile type.
    kperk_pairs : iterable
        an iterable of pairs of kper (zero-based stress
        period index) and k (zero-based layer index) to
        setup observations for.  If None, then a shit-ton
        of observations may be produced!
    skip : variable
        a value or function used to determine which values
        to skip when setting up observations.  If np.scalar(skip)
        is True, then values equal to skip will not be used.
        If skip can also be a np.ndarry with dimensions equal to the model.
        Obscervations are set up only for cells with Non-zero values in the array.
        If not np.ndarray or np.scalar(skip), then skip will be treated as a lambda function that
        returns np.NaN if the value should be skipped.
    prefix : str
        the prefix to use for the observation names. default is "hds".

    Returns
    -------
    (forward_run_line, df) : str, pd.DataFrame
        a python code str to add to the forward run script and the setup info for the observations

    Note
    ----
    requires flopy

    writes <hds_file>.dat.ins instruction file

    writes _setup_<hds_file>.csv which contains much
    useful information for construction a control file


    """
    try:
        import flopy
    except Exception as e:
        print("error importing flopy, returning {0}".format(str(e)))
        return

    assert os.path.exists(hds_file),"head save file not found"
    if hds_file.lower().endswith(".ucn"):
        try:
            hds = flopy.utils.UcnFile(hds_file)
        except Exception as e:
            raise Exception("error instantiating UcnFile:{0}".format(str(e)))
    else:
        try:
            hds = flopy.utils.HeadFile(hds_file)
        except Exception as e:
            raise Exception("error instantiating HeadFile:{0}".format(str(e)))

    if kperk_pairs is None:
        kperk_pairs = []
        for kstp,kper in hds.kstpkper:
            kperk_pairs.extend([(kper-1,k) for k in range(hds.nlay)])
    if len(kperk_pairs) == 2:
        try:
            if len(kperk_pairs[0]) == 2:
                pass
        except:
            kperk_pairs = [kperk_pairs]

    #if start_datetime is not None:
    #    start_datetime = pd.to_datetime(start_datetime)
    #    dts = start_datetime + pd.to_timedelta(hds.times,unit='d')
    data = {}
    kpers = [kper-1 for kstp,kper in hds.kstpkper]
    for kperk_pair in kperk_pairs:
        kper,k = kperk_pair
        assert kper in kpers, "kper not in hds:{0}".format(kper)
        assert k in range(hds.nlay), "k not in hds:{0}".format(k)
        kstp = last_kstp_from_kper(hds,kper)
        d = hds.get_data(kstpkper=(kstp,kper))[k,:,:]

        data["{0}_{1}".format(kper,k)] = d.flatten()
        #data[(kper,k)] = d.flatten()
    idx,iidx,jidx = [],[],[]
    for _ in range(len(data)):
        for i in range(hds.nrow):
            iidx.extend([i for _ in range(hds.ncol)])
            jidx.extend([j for j in range(hds.ncol)])
            idx.extend(["i{0:04d}_j{1:04d}".format(i,j) for j in range(hds.ncol)])
    idx = idx[:hds.nrow*hds.ncol]

    df = pd.DataFrame(data,index=idx)
    data_cols = list(df.columns)
    data_cols.sort()
    #df.loc[:,"iidx"] = iidx
    #df.loc[:,"jidx"] = jidx
    if skip is not None:
        for col in data_cols:
            if np.isscalar(skip):
                df.loc[df.loc[:,col]==skip,col] = np.NaN
            elif isinstance(skip, np.ndarray):
                assert skip.ndim >= 2, "skip passed as {}D array, At least 2D (<= 4D) array required".format(skip.ndim)
                assert skip.shape[-2:] == (hds.nrow, hds.ncol), \
                    "Array dimensions of arg. skip needs to match model dimensions ({0},{1}). ({2},{3}) passed".\
                        format(hds.nrow, hds.ncol, skip.shape[-2], skip.shape[-1])
                if skip.ndim == 2:
                    print("2D array passed for skip, assuming constant for all layers and kper")
                    skip = np.tile(skip, (len(kpers), hds.nlay, 1, 1))
                if skip.ndim == 3:
                    print("3D array passed for skip, assuming constant for all kper")
                    skip = np.tile(skip, (len(kpers), 1, 1, 1))
                kper, k = [int(c) for c in col.split('_')]
                df.loc[df.index.map(
                    lambda x: skip[kper, k, int(x.split('_')[0].strip('i')), int(x.split('_')[1].strip('j'))] == 0),
                    col] = np.NaN
            else:
                df.loc[:,col] = df.loc[:,col].apply(skip)

    # melt to long form
    df = df.melt(var_name="kperk",value_name="obsval")
    # set row and col identifies
    df.loc[:,"iidx"] = iidx
    df.loc[:,"jidx"] = jidx
    #drop nans from skip
    df = df.dropna()
    #set some additional identifiers
    df.loc[:,"kper"] = df.kperk.apply(lambda x: int(x.split('_')[0]))
    df.loc[:,"kidx"] = df.pop("kperk").apply(lambda x: int(x.split('_')[1]))

    # form obs names
    #def get_kper_str(kper):
    #    if start_datetime is not None:
    #        return  dts[int(kper)].strftime("%Y%m%d")
    #    else:
    #        return "kper{0:04.0f}".format(kper)
    fmt = prefix + "_{0:02.0f}_{1:03.0f}_{2:03.0f}_{3:03.0f}"
    # df.loc[:,"obsnme"] = df.apply(lambda x: fmt.format(x.kidx,x.iidx,x.jidx,
    #                                                   get_kper_str(x.kper)),axis=1)
    df.loc[:,"obsnme"] = df.apply(lambda x: fmt.format(x.kidx,x.iidx,x.jidx,
                                                      x.kper),axis=1)

    df.loc[:,"ins_str"] = df.obsnme.apply(lambda x: "l1 w !{0}!".format(x))
    df.loc[:,"obgnme"] = prefix
    #write the instruction file
    with open(hds_file+".dat.ins","w") as f:
        f.write("pif ~\nl1\n")
        df.ins_str.to_string(f,index=False,header=False)

    #write the corresponding output file
    df.loc[:,["obsnme","obsval"]].to_csv(hds_file+".dat",sep=' ',index=False)

    hds_path = os.path.dirname(hds_file)
    setup_file = os.path.join(hds_path,"_setup_{0}.csv".format(os.path.split(hds_file)[-1]))
    df.to_csv(setup_file)
    fwd_run_line = "pyemu.gw_utils.apply_hds_obs('{0}')\n".format(hds_file)
    df.index = df.obsnme
    return fwd_run_line, df