diff --git a/bin/schemas/snaphu_schema.yaml b/bin/schemas/snaphu_schema.yaml
new file mode 100644
index 0000000..cb0c8b9
--- /dev/null
+++ b/bin/schemas/snaphu_schema.yaml
@@ -0,0 +1,409 @@
+snaphu:
+    # REQUIRED. Path to unwrapped phase rasters. GDT_Float32,same shape as igram_raster.
+    unwrapped_raster: list(str(), required=False)
+
+    # REQUIRED. Path to connected components rasters. GDT_UInt32, same shape as igram_raster
+    connected_components_raster: list(str(), required=False)
+
+    # REQUIRED. Path to input wrapped interferograms (GDT_CFloat32)
+    igram_raster: list(str())
+
+    # REQUIRED. Path to correlation magnitudes normalized in [0, 1]. GDT_Float32, same shape as igram_raster
+    correlation_raster: list(str())
+
+    # REQUIRED Effective number of looks used in input correlation data
+    nlooks: num(min=1)
+
+    # Statistical cost mode
+    cost_mode: enum('topo', 'defo', 'smooth', 'p-norm', required=False)
+
+    # Configuration parameters for the specified cost mode.
+    cost_mode_parameters: include('cost_mode_options', required=False)
+
+    # File path to binary mask of valid pixels (GDT_Byte). Pixels labelled
+    # with 0 will be masked out in the wrapped interferogram.
+    mask: list(str(), required=False)
+
+    # Coarse unwrapped phase estimates to guide SNAPHU. Must be same shape
+    # as input wrapped interferogram.
+    unwrapped_phase_estimate: list(str(), required=False)
+
+    # Initialization method (MST: Minimum Spanning tree;
+    # MCF: Minimum Cost Flow)
+    initialization_method: enum('mst', 'mcf', required=False)
+
+    # Average intensity of the two SLCs in linear units (not dB). If None, interfeogram magnitude
+    # is used as intensity. Used only in topo cost mode!
+    power: list(str(), required=False)
+
+    # Parameters for SNAPHU tiling mode
+    tiling_parameters: include('tiling_options', required=False)
+
+    # Configuration parameters used by the network initialization and nonlinear
+    # network flow solver algorithms
+    solver_parameters: include('solver_options', required=False)
+
+    # Connected components parameters
+    connected_components_parameters: include('connected_components_options', required=False)
+
+    # Parameters for correlation bias model
+    correlation_bias_parameters: include('correlation_bias_options', required=False)
+
+    # Parameters for phase standard deviation model
+    phase_std_parameters: include('phase_std_options', required=False)
+
+    # Path to directory to store intermediate results. If not specified or
+    # not exists,a temporary directory is created and automatically removed to
+    # from the file system at the end of the processing. Otherwise, the directory
+    # and its content will not be cleaned up
+    scratchdir: str(required=False)
+
+    # If True,prints verbose outputs (default: False)
+    verbose: bool(required=False)
+
+    # If True, dumps intermediate data to 'scratchdir' for debugging
+    debug: bool(required=False)
+
+---
+cost_mode_options:
+    # Topo cost mode parameters
+    topo_parameters: include('topo_cost_options', required=False)
+
+    # Deformation cost parameter options
+    deformation_parameters: include('deformation_cost_options', required=False)
+
+    # Smooth cost parameter options
+    smooth_parameters: include('smooth_cost_options', required=False)
+
+    # PNorm cost parameter options
+    pnorm_parameters: include('pnorm_cost_options', required=False)
+
+
+topo_cost_options:
+    # Note: bperp, near_range, dr, da, range_res, azimuth_res, wavelength,
+    # transmit_mode, and altitude are set as required because needed by
+    # TopoCostParams constructor
+
+    # Perpendicular baseline length (meters). Negative values imply
+    # increasing topographic heights.
+    bperp: num()
+
+    # Slant range from platform to first range bin (meters)
+    near_range: num()
+
+    # Slant range bin spacing after multi-looking (meters)
+    dr: num(min=0)
+
+    # Azimuth bin spacing after multi-looking (meters)
+    da: num(min=0)
+
+    # Single-look slant range resolution (meters)
+    range_res: num(min=0)
+
+    # Single-look azimuth resolution (meters)
+    azimuth_res: num(min=0)
+
+    # Wavelength (meters)
+    wavelength: num(min=0)
+
+    # Radar transmit mode. For 'pingpong' and 'repeat_pass' indicate that both
+    # antennas transmitted and received (same effect on the algorithm).
+    # In 'single_antenna_transmit' mode a single antenna was used to transmit
+    # while both antennas received. In this mode, the baseline is effectively halved.
+    transmit_mode: enum('pingpong', 'repeat_pass', 'single_antenna_transmit')
+
+    # Sensor platform altitude relative to the Earth's surface (meters)
+    altitude: num(min=0)
+
+    # Local Earth's radius (meters) for spherical Earth model (default: 6378000.0)
+    earth_radius: num(min=0, required=False)
+
+    # Ratio of diffuse to specular scattering (default: 0.02)
+    kds: num(required=False)
+
+    # Power specular scattering component. Large values imply larger peak
+    # for specular scattering (default: 8.0)
+    specular_exp: num(required=False)
+
+    # Multiplicative factor applied to diffuse scatter term to evaluate crossover
+    # point between diffuse and specular scatter in terms of range slope (default: 2.0)
+    dzr_crit_factor: num(required=False)
+
+    # Flag to enable discontinuities from shadowing. If False, the min topographic
+    # slope is estimated from mean backscatter intensity clipped to 'dz_ei_min'. (default: False)
+    shadow: bool(required=False)
+
+    # Min slope expected in absence of layover (meters per slant range pixel, default: -4.0)
+    dz_ei_min: num(required=False)
+
+    # Width of window (number of pixels) for summing layover brightness (default: 16)
+    layover_width: int(min=0, required=False)
+
+    # Normalized threshold brightness for assuming layover (default: 1.25)
+    layover_min_ei: num(required=False)
+
+    # Multiplicative factor applied to kds to get ratio of slopes for linearized
+    # scattering model. The term improves agreement of piecewise-linear model with
+    # the cosine model near the transition point (dzrcrit) at the expense of poorer
+    # agreement at very large slopes (default: 1.18)
+    slope_ratio_factor: num(required=False)
+
+    # Variance of range slopes due to uncertainties in slope estimation from
+    # brightness in (meter/pixels)^2 (default: 100.0)
+    sigsq_ei: num(required=False)
+
+    # Step size for lookup table of max layover slope based on measured
+    # normalized correlation (default: 0.005)
+    drho: num(required=False)
+
+    # Layover peak location (meters/pixels, default: -2.0)
+    dz_layover_peak: num(required=False)
+
+    # Factor applied to range layover probability to get azimuth layover
+    # probability density (default: 0.99)
+    azdz_factor: num(required=False)
+
+    # Factor applied to slope expected from brightness without layover. Can
+    # account for underestimation of brightness from averaging with neighboring
+    # dark pixels when despeckling. (default: 4.0)
+    dz_ei_factor: num(required=False)
+
+    # eight applied to slope expected from brightness without layover. Must
+    # be between zero and one. Can reduce influence of intensity on
+    # on-layover slope. This is useful if there are lots of non-topographic
+    # variations in brightness (i.e. changes in surface reflectivity, default: 0.5)
+    dz_ei_weight: num(required=False)
+
+    # Factor applied to slope expected from brightness with layover. Can
+    # account for underestimation of brightness from averaging with
+    # neighboring dark pixels when despeckling. (default: 1.0)
+    dz_layover_factor: num(required=False)
+
+    # Ratio of layover probability density to peak probability density for
+    # non-layover slopes expected. (default: 0.9)
+    layover_const: num(required=False)
+
+    # Factor applied to slope variance for non-layover to get falloff of
+    # probability density after the  upper layover slope limit has been
+    # exceeded. (default: 2.0)
+    layover_falloff_const: num(required=False)
+
+    # Fraction of (ambiguity height)^2 to use for slope variance in the
+    # presence of layover. (default: 0.1)
+    sigsq_layover_factor: num(required=False)
+
+    # Number of rows to use in sliding average window for normalized
+    # intensity values (default: 65)
+    krow_ei: int(min=0, required=False)
+
+    # Same as above but for columns (default: 257)
+    kcol_ei: int(min=0, required=False)
+
+    # Initial value of range slope for dzrcrit numerical solution
+    # (meters/pixel, default: 2048.0)
+    init_dzr: num(required=False)
+
+    # Initial range slope step size in dzrhomaz numerical solution
+    # (meters/pixel, default: 100.0)
+    init_dz_step: num(required=False)
+
+    # Height of ambiguity for auto-scaling the `SolverParams.cost_scale`
+    # parameter to equal 100. This is the amount of height change, in meters,
+    # resulting in a :math:`2 \pi` change in the interferometric phase. The
+    # cost scale is automatically adjusted to be inversely proportional to the
+    # mid-swath ambiguity height. (default: 80.0)
+    cost_scale_ambiguity_height: num(required=False)
+
+    # Step size, in radians, for dzrhomax lookup table. The index is on the
+    # flat-earth incidence angle; this is the sample spacing in the table.
+    # (default: 0.01)
+    dnom_inc_angle: num(required=False)
+
+    # Number sliding window pixel to average the wrapped phase gradient to get mean
+    # non-layover slope in the direction parallel to the examined phase difference (default: 7)
+    kpar_dpsi: int(min=1, required=False)
+
+    # Same as above but in the perpendicular direction
+    kperp_dpsi: int(min=1, required=False)
+
+deformation_cost_options:
+    # Factor applied to range discontinuity probability density to get
+    # corresponding value for azimuth (default: 1.0)
+    azdz_factor: num(required=False)
+
+    # Max phase discontinuity in units of cycle of 2*pi. Set to 0 if not abrupt
+    # phase discontinuity are expected (default: 1.2)
+    defo_max: num(min=0, required=False)
+
+    # Phase variance in cycle^2 reflecting uncertainty in measurement of
+    # actual statistical correlation (default: 0.05)
+    sigsq_corr: num(min=0, required=False)
+
+    # Ratio of phase discontinuity probability density to peak probability
+    # density expected for discontinuity-possible pixel differences. 1: zero
+    # cost for discontinuity, o: infinite cost (default: 0.9)
+    defo_const: num(min=0, max=1, required=False)
+
+    # Factor for slope variance for non-layover to get falloff of probability
+    # density after the upper layover slope limit has been exceeded (default: 2.0)
+    layover_falloff_const: num(required=False)
+
+    # Number sliding window pixel to average the wrapped phase gradient to get mean
+    # non-layover slope in the direction parallel to the examined phase difference (default: 7)
+    kpar_dpsi: int(min=1, required=False)
+
+    # Same as above but in the perpendicular direction (default: 7)
+    kperp_dpsi: int(min=1, required=False)
+
+smooth_cost_options:
+    # Number sliding window pixel to average the wrapped phase gradient to get mean
+    # non-layover slope in the direction parallel to the examined phase difference (default: 7)
+    kpar_dpsi: int(min=1, required=False)
+
+    # Same as above but in the perpendicular direction (default: 7)
+    kperp_dpsi: int(min=1, required=False)
+
+pnorm_cost_options:
+    # Lp norm non-negative exponent.
+    lp_exp: num(min=0, required=False)
+
+    # Flag to activate bidirectional Lp costs. If True, scalar weight of a Lp arc
+    # may be different depending on the direction of net flow arc. If False, scalar
+    # weight is the same regardless arc direction
+    bidirection: bool(required=False)
+
+tiling_options:
+    # Max number of child processes to spawn for parallel tile unwrapping. If < 1
+    # use all available processors (default: 1)
+    nproc: int(required=False)
+
+    # Number of tiles along row direction. If `tile_nrows`=`tile_ncols` = 1,
+    # the interferogram is unwrapped as a single (default: 1)
+    tile_nrows: int(min=1, required=False)
+
+    # Same as above but for columns (default: 1)
+    tile_ncols: int(min=1, required=False)
+
+    # Number of overlapping rows between neighboring tiles (default: 0)
+    row_overlap: int(min=0, required=False)
+
+    # Number of overlapping columns between neighboring tiles (default: 0)
+    col_overlap: int(min=0, required=False)
+
+    # Cost threshold for determining boundaries of reliable regions.
+    # Larger cost threshold implies smaller regions (safer, but more expensive
+    # computationally). (default: 500)
+    tile_cost_thresh: int(required=False)
+
+    # Min size (in pixels) of a reliable region in tile mode(default: 100)
+    min_region_size: int(required=False)
+
+    # Extra weight applied to secondary arcs on tile edges (default: 2.5)
+    tile_edge_weight: num(required=False)
+
+    # Max flow magnitude whose cost will be stored in the secondary cost LUT.
+    # Secondary costs larger than this will be approximated by a quadratic function. (default: 8)
+    secondary_arc_flow_max: int(required=False)
+
+    # If True, re-optimize as a single tile after using tile mode for
+    # initialization. Equivalent to unwrapping with multiple tiles,
+    # then using the unwrapped output as the input to a new, SNAPHU single-tile run
+    # making iterative improvements to the solution. May improve speed compared to
+    # a single single-tile run. (default: False)
+    single_tile_reoptimize: bool(required=False)
+
+solver_options:
+    # Maximum flow increment (default: 4)
+    max_flow_increment: int(min=0, required=False)
+
+    # Maximum flow in initialization. If 0, max flow is computed automatically.
+    # To disable, set to large value not overflowing long integer data (default: 9999)
+    initial_max_flow: int(min=0, max=9999, required=False)
+
+    # Constant to add to max flow expected from statistical cost functions for automatically
+    # determining initial maximum flow (default: 3)
+    arc_max_flow_const: int(required=False)
+
+    # Threshold precision for iterative numerical calculation (default: 0.001)
+    threshold: num(required=False)
+
+    # Max cost for scalar MST cost and to estimate number of buckets
+    # for solver routine (default: 1000.0)
+    max_cost: num(required=False)
+
+    # Cost scaling factor applied to floating-point costs before integer
+    # cost quantization (default: 100.0)
+    cost_scale: num(required=False)
+
+    # Integer spacing representing one unit of flow (one phase cycle)
+    # when storing costs as short integers (default: 200)
+    n_cycle: int(required=False)
+
+    # Fraction of total number of nodes to add in each tree expansion phase of
+    # the solver algorithm. (default: 0.0008)
+    max_new_node_const: num(required=False)
+
+    # Number of cycles to call to the solver with a specific flow increment delta
+    # and still consider that increment done. Ideally zero, but scaling for different
+    # deltas may leave some negative cycles not critically affecting the solution.
+    # If None, automatically determined from interferogram shape (default: None)
+    max_n_flow_cycles: num(required=False)
+
+    # Fraction of the number of pixels to use as the maximum number of cycles
+    # allowed for a specific flow increment if max_n_flow_cycles was None (default: 0.00001)
+    max_cycle_frac: num(required=False)
+
+    # Minimum number of connected nodes for unwrappping. If disconnected sets of pixels are
+    # present, a source is selected for each connected set, if number of nodes in the set
+    # is greater than n_conn_node_min.
+    n_conn_node_min: int(min=0, required=False)
+
+    # Number of major iteration between tree pruning operations. Smaller numbers cause
+    # pruning to occur more frequently (default: 2000000000)
+    n_major_prune: int(required=False)
+
+    # Cost threshold for tree pruning. Lower thresholds prune more aggressively
+    # (default: 2000000000)
+    prune_cost_threshold: int(required=False)
+
+connected_components_options:
+    # Minimum size of a single connected component, as a fraction of the total
+    # number of pixels in the tile (default: 0.01)
+    min_frac_area: num(required=False)
+
+    # Cost threshold for connected components. Higher threshold gives smaller
+    # connected components (default: 300)
+    cost_threshold: int(required=False)
+
+    # Max number of connected components per tile (default: 32)
+    max_ncomps: int(required=False)
+
+correlation_bias_options:
+    # c1 parameter for correlation bias model (Touzi et al. 1999)
+    c1: num(required=False)
+
+    # c2 parameter for correlation bias model (Touzi et al. 1999)
+    c2: num(required=False)
+
+    # Factor applied to expected minimum measured (biased) correlation
+    # coefficient. Values smaller than the threshold min_corr_factor * rho0
+    # are assumed to come from zero statistical correlation because of
+    # estimator bias. rho0 is the expected biased correlation measure if the
+    # true correlation is zero. (default: 1.25)
+    min_corr_factor: num(required=False)
+
+# Note: Interferometric phase standard deviation is modelled as
+# .. math:: \sigma_{\phi} = \rho ^{ c_1 + c_2 * \log nlooks + c_3 * nlooks }
+phase_std_options:
+    # c1 parameter for phase std model
+    c1: num(required=False)
+
+    # c2 parameter for phase std model
+    c2: num(required=False)
+
+    # c3 parameter for phase std model
+    c3: num(required=False)
+
+    # Min phase variance value after quantization to integers.
+    # Must be greater than 0 to prevent division by 0. (default: 1)
+    sigsq_min: int(required=False)
\ No newline at end of file
diff --git a/bin/unwrap_snaphu.py b/bin/unwrap_snaphu.py
new file mode 100644
index 0000000..782bad8
--- /dev/null
+++ b/bin/unwrap_snaphu.py
@@ -0,0 +1,700 @@
+import argparse
+import os
+import isce3
+import isce3.unwrap.snaphu as snaphu
+import yamale
+import journal
+import numpy as np
+from ruamel.yaml import YAML
+
+EXAMPLE = """example:
+  unwrap_snaphu.py <custom_runconfig>   # Run with default and custom runconfigs
+  unwrap_snaphu.py -h / --help          # help
+"""
+
+WORKFLOW_SCRIPTS_DIR = os.path.dirname(os.path.realpath(__file__))
+
+def command_line_parser():
+    """
+    Command line parser
+    """
+    parser = argparse.ArgumentParser(
+        description="Unwrap interferograms using SNAPHU unwrapper",
+        formatter_class=argparse.RawTextHelpFormatter,
+        epilog=EXAMPLE)
+    parser.add_argument('-r', '--custom_runconfig', type=str,
+                        dest='run_config_path',
+                        help='YAML custom runconfig for SNAPHU unwrapping. '
+                             'Ignored if default runconfig snaphu.yaml is input.')
+    return parser.parse_args()
+
+
+def load_runconfig(cmd_opts):
+    """
+    Load SNAPHU runconfig. If cmd_opts.run_config_path
+    is None, loads default runconfig. Otherwise, default
+    runconfig is updated with user-defined options.
+    Parameters
+    ----------
+    cmd_opts: NameSpace
+       Command line option parser
+    """
+    error_channel = journal.error('unwrap_snaphu.load_runconfig')
+    try:
+        # Load schemas corresponding to SNAPHU and to validate against
+        schema = yamale.make_schema(f'{WORKFLOW_SCRIPTS_DIR}/schemas/snaphu_schema.yaml',
+                                    parser='ruamel')
+    except:
+        err_str = f'Unable to load the schema for SNAPHU unwrapper'
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
+    runconfig_path = cmd_opts.run_config_path
+    if os.path.isfile(runconfig_path):
+        try:
+            data = yamale.make_data(runconfig_path, parser='ruamel')
+        except yamale.YamaleError as yamale_err:
+            err_str = f'Yamale is unable to load the SNAPHU runconfig at {runconfig_path}'
+            error_channel.log(err_str)
+            raise yamale.YamaleError(err_str) from yamale_err
+    else:
+        err_str= f'Runconfig file at {runconfig_path} has not been found'
+        error_channel.log(err_str)
+        raise FileNotFoundError(err_str)
+
+    # Validate YAML file from command line
+    try:
+        yamale.validate(schema, data)
+    except yamale.YamaleError as yamale_err:
+        err_str = f'Schema validation failed for SNAPHU and runconfig at {runconfig_path}'
+        error_channel.log(err_str)
+        raise yamale.YamaleError(err_str) from yamale_err
+
+    # Load user-provided runconfig
+    parser = YAML(typ='safe')
+    with open(runconfig_path, 'r') as f:
+        runconfig = parser.load(f)
+    return runconfig
+
+
+def clean_none_cfg_keys(cfg):
+    """
+    Clean dict from None values
+    Parameters
+    ----------
+    cfg: dict
+      Dictionary to be cleaned
+    Returns
+    -------
+    clean_cfg: dict
+      Dictionary without None values
+    """
+    clean_cfg = {}
+    for k, v in cfg.items():
+        if isinstance(v, dict):
+            nested = clean_none_cfg_keys(v)
+            if len(nested.keys()) > 0:
+                clean_cfg[k] = nested
+        elif v is not None:
+            clean_cfg[k] = v
+    return clean_cfg
+
+
+def check_dataset_group(cfg, igram_len, name):
+    """
+    Check ancillary raster count. Ancillary rasters
+    are mask, power, unwrapped_phase_estimate. User
+    can provide one raster or as many raster as the
+    number of interferograms to unwrapped.
+    cfg: dict
+       Dictionary containing SNAPHU paramters
+    igram_len: int
+       Number of interferograms to unwrap
+    name: str
+       Name of dataset to check (e.g. 'mask')
+    """
+    error_channel = journal.error(
+        'main.check_input_files.check_dataset_group')
+
+    if (cfg.get(name) is not None):
+        data_cfg = cfg[name]
+        if len(data_cfg) > 1:
+           if len(data_cfg) != igram_len:
+               err_str = f'Assign one {name} or as many {name}s as number of wrapped igrams'
+               error_channel.log(err_str)
+           for file in data_cfg:
+               if not os.path.exists(file):
+                   err_str = f'{file} for dataset {name} does not exist'
+                   error_channel.log(err_str)
+
+
+def check_input_files(cfg):
+    """
+    Check input file group of SNAPHU runconfig
+    Parameters
+    ----------
+    cfg: dict
+      Dictionary containing SNAPHU parameters
+    """
+    error_channel = journal.error('main.check_input_files')
+    # Check that input cfg is valid (not None)
+    if cfg is None:
+        err_str = 'SNAPHU parameter dictionary is not valid (None)'
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
+    # Extract input & output files to check
+    igram_paths = cfg['igram_raster']
+    corr_paths = cfg['correlation_raster']
+    ugram_paths = cfg['unwrapped_raster']
+    conn_comp_paths = cfg['connected_components_raster']
+
+    # Check consistency between number of interferograms and correlations
+    if len(igram_paths) != len(corr_paths):
+        err_str = "Number of interferograms must be consistent with number of correlations"
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
+    # Check that user-defined interferograms and correlations exist
+    for file_igram, file_corr in zip(igram_paths, corr_paths):
+        if not os.path.exists(file_igram) or not os.path.exists(file_corr):
+            err_str = f'{file_igram} or {file_corr} files do not exists'
+            error_channel.log(err_str)
+            raise FileNotFoundError(err_str)
+
+    # If assigned, number of unwrapped/connected components file paths
+    # must correspond to number of interferogram file paths
+    if (ugram_paths is not None) and (len(ugram_paths) != len(igram_paths)):
+        err_str = f'Number of unwrapped file paths must correspond to number of igram file paths'
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+    if (conn_comp_paths is not None) and (
+            len(conn_comp_paths) != len(igram_paths)):
+        err_str = f'Number of connected components file paths must correspond to number of igram file paths'
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
+    # For mask, power, and unwrapped phase estimate, use one raster for
+    # all the interferograms or as many rasters as the number of interferograms
+    check_dataset_group(cfg, len(igram_paths), 'mask')
+    check_dataset_group(cfg, len(igram_paths), 'unwrapped_phase_estimate')
+    check_dataset_group(cfg, len(igram_paths), 'power')
+
+
+def set_topo_params(cost_cfg):
+    """
+    Set parameters for SNAPHU topo cost mode
+    Parameters
+    ----------
+    cost_cfg: dict
+      Dictionary containing SNAPHU cost parameters options
+    Returns
+    -------
+    topo: isce3.unwrap.snaphu.TopoCostParams
+       Topo cost parameter object with user-defined parameters
+    """
+    # Note: isce3.snaphu.unwrap.TopoCostParams is a frozen dataclass.
+    # Therefore, attributes are not settable. We use dict.get to initialize object
+    # dict.get() will assign the value stored in dict keyword if keyword is found
+    # or a specified default value. More on assigned default values at:
+    # https://github-fn.jpl.nasa.gov/isce-3/isce/blob/develop/python/packages/isce3/unwrap/snaphu.py
+
+    # Note: bperp, near_range, dr, dz, range_res, az_red, wavelength, transmit_mode
+    # altitude are required to initialize object. Schema will check their existence
+
+    error_channel = journal.error('unwrap_snaphu.set_topo_params')
+
+    if 'topo_parameters' in cost_cfg:
+        cfg = cost_cfg['topo_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        topo = snaphu.TopoCostParams(bperp=cfg.get('bperp'),
+                                     near_range=cfg.get('near_range'),
+                                     dr=cfg.get('dr'), da=cfg.get('da'),
+                                     range_res=cfg.get('range_res'),
+                                     az_res=cfg.get('azimuth_res'),
+                                     wavelength=cfg.get('wavelength'),
+                                     transmit_mode=cfg.get('transmit_mode'),
+                                     altitude=cfg.get('altitude'),
+                                     earth_radius=cfg.get('earth_radius',
+                                                          6378000.0),
+                                     kds=cfg.get('kds', 0.02),
+                                     specular_exp=cfg.get('specular_exp', 8.0),
+                                     dzr_crit_factor=cfg.get('dzr_crit_factor',
+                                                             2.0),
+                                     shadow=cfg.get('shadow', False),
+                                     dz_ei_min=cfg.get('dz_ei_min', -4.0),
+                                     lay_width=cfg.get('layover_width', 16),
+                                     lay_min_ei=cfg.get('layover_min_ei', 1.25),
+                                     slope_ratio_factor=cfg.get(
+                                         'slope_ratio_factor', 1.18),
+                                     sigsq_ei=cfg.get('sigsq_ei', 100.0),
+                                     drho=cfg.get('drho', 0.005),
+                                     dz_lay_peak=cfg.get('dz_layover_peak',
+                                                         -2.0),
+                                     azdz_factor=cfg.get('azdz_factor', 0.99),
+                                     dz_ei_factor=cfg.get('dz_ei_factor', 4.0),
+                                     dz_ei_weight=cfg.get('dz_ei_weight', 0.5),
+                                     dz_lay_factor=cfg.get('dz_layover_factor',
+                                                           1.0),
+                                     lay_const=cfg.get('layover_const', 0.9),
+                                     lay_falloff_const=cfg.get(
+                                         'layover_falloff_const', 2.0),
+                                     sigsq_lay_factor=cfg.get(
+                                         'sigsq_layover_factor', 0.1),
+                                     krow_ei=cfg.get('krow_ei', 65),
+                                     kcol_ei=cfg.get('kcol_ei', 257),
+                                     init_dzr=cfg.get('init_dzr', 2048.0),
+                                     init_dz_step=cfg.get('init_dz_step', 100),
+                                     cost_scale_ambig_ht=cfg.get(
+                                         'cost_scale_ambiguity_height', 80.0),
+                                     dnom_inc_angle=cfg.get('dnom_inc_angle',
+                                                            0.01),
+                                     kpar_dpsi=cfg.get('kpar_dpsi', 7),
+                                     kperp_dpsi=cfg.get('kperp_dpsi', 7))
+    else:
+        err_str = "bperp, near_range, dr, dz, range_res, az_red, wavelength, " \
+                  "transmit_mode, altitude are required to initialize topo cost parameter object "
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+    return topo
+
+
+def set_defo_params(cost_cfg):
+    """
+    Set parameters for SNAPHU deformation cost mode
+    Parameters
+    ----------
+    cost_cfg: dict
+      Dictionary containing SNAPHU cost parameters
+    Returns
+    -------
+    defo: isce3.unwrap.snaphu.DefoCostParams
+      Deformation cost parameter object with user-defined parameters
+    """
+    # If deformation parameter is not empty, extract it
+    if cost_cfg.get('deformation_parameters') is not None:
+        cfg = cost_cfg['deformation_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        defo = snaphu.DefoCostParams(azdz_factor=cfg.get('azdz_factor', 1.0),
+                                     defo_max=cfg.get('defo_max', 1.2),
+                                     sigsq_corr=cfg.get('sigsq_corr', 0.05),
+                                     defo_const=cfg.get('defo_const', 0.9),
+                                     lay_falloff_const=cfg.get(
+                                         'layover_falloff_const', 2.0),
+                                     kpar_dpsi=cfg.get('kpar_dpsi', 7),
+                                     kperp_dpsi=cfg.get('kperp_dpsi', 7))
+    else:
+        defo = snaphu.DefoCostParams()
+    return defo
+
+
+def set_smooth_params(cost_cfg):
+    """
+    Set parameters for SNAPHU smooth cost mode
+    Parameters
+    ----------
+    cost_cfg: dict
+      Dictionary containing SNAPHU cost parametersoptions
+    Returns
+    -------
+    smooth: isce3.unwrap.snaphu.SmoothCostParams
+      Smooth cost parameter object with user-defined parameters
+    """
+    # If smooth parameters are present, extract smooth dict
+    if cost_cfg.get('smooth_parameters') is not None:
+        cfg = cost_cfg['smooth_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        smooth = snaphu.SmoothCostParams(kpar_dpsi=cfg.get('kpar_dpsi', 7),
+                                         kperp_dpsi=cfg.get('kperp_dpsi', 7))
+    else:
+        # use all defaults values
+        smooth = snaphu.SmoothCostParams()
+    return smooth
+
+
+def set_pnorm_params(cost_cfg):
+    """
+    Set parameters for SNAPHU P-Norm cost mode
+    Parameters
+    ----------
+    cost_cfg: dict
+      Dictionary containing SNAPHU cost parameter options
+    Returns
+    -------
+    pnorm: isce3.unwrap.snaphu.PNormCostParams
+      P-Norm cost parameter object with user-defined parameters
+    """
+
+    # If pnorm section of runconfig is not empty,
+    # proceed to set user-defined pnorm parameters
+    if cost_cfg.get('pnorm_parameters') is not None:
+        cfg = cost_cfg['pnorm_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        pnorm = snaphu.PNormCostParams(p=cfg.get('lp_exp', 0.0),
+                                       bidir=cfg.get('bidirection', True))
+    else:
+        pnorm = snaphu.PNormCostParams()
+    return pnorm
+
+
+def select_cost_options(cfg, cost_mode='defo'):
+    """
+    Select and set cost parameter object based
+    on user-defined cost mode options
+    Parameters
+    ----------
+    cfg: dict
+      Dictionary containing SNAPHU parameter options
+    cost_mode: str
+      Snaphu cost mode. Default: defo
+    Returns
+    -------
+    cost_params_obj: object
+      Object corresponding to selected cost mode.
+      E.g. if cost_mode is "defo", cost_params_obj
+      is an instance of isce3.unwrap.snaphu.DefoCostParams()
+    """
+    error_channel = journal.error('unwrap_snaphu.select_cost_options')
+
+    # If 'cost_mode_parameters' does not exist, create empty dictionary
+    if cfg.get('cost_mode_parameters') is None:
+        cfg['cost_mode_parameters'] = {}
+
+    if cost_mode == 'topo':
+        cost_params_obj = set_topo_params(cfg['cost_mode_parameters'])
+    elif cost_mode == 'defo':
+        cost_params_obj = set_defo_params(cfg['cost_mode_parameters'])
+    elif cost_mode == 'smooth':
+        cost_params_obj = set_smooth_params(cfg['cost_mode_parameters'])
+    elif cost_mode == 'p-norm':
+        cost_params_obj = set_pnorm_params(cfg['cost_mode_parameters'])
+    else:
+        err_str = f"{cost_mode} is not a valid cost mode option"
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+    return cost_params_obj
+
+
+def get_raster(cfg, name, iter=0):
+    """
+    Select ancillary raster (e.g. mask) based
+    on dataset name and number.
+    Parameters:
+    ----------
+    cfg: dict
+      Dictionary containing SNAPHU parameter options
+    name: str
+      Name of the ancillary raster to extract
+    iter: int
+      Iterator to extract correct dataset
+    """
+    if cfg.get(name) is not None:
+        # If dataset exist, check number of allocated rasters
+        if len(cfg[name]) == 1:
+            path = cfg[name]
+        # If more than one, use iter to extract correct raster
+        else:
+            path = cfg[name][iter]
+        raster = isce3.io.gdal.Raster(path)
+    else:
+        raster = None
+    return raster
+
+
+def set_tile_params(snaphu_cfg):
+    """
+    Set user-defined tiling parameter options
+    Parameters
+    ----------
+    snaphu_cfg: dict
+       Dictionary containing SNAPHU options
+    Returns:
+    -------
+    tile: isce3.unwrap.snaphu.TilingParams() or None
+       Object containing tiling parameters options or None
+    """
+    # If 'tiling_parameters' is in snaphu_cfg, inspect setted
+    # options. If None found, assign default
+    if snaphu_cfg.get('tiling_parameters') is not None:
+        cfg = snaphu_cfg['tiling_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        tile = snaphu.TilingParams(nproc=cfg.get('nproc', 1),
+                                   tile_nrows=cfg.get('tile_nrows', 1),
+                                   tile_ncols=cfg.get('tile_ncols', 1),
+                                   row_overlap=cfg.get('row_overlap', 0),
+                                   col_overlap=cfg.get('col_overlap', 0),
+                                   tile_cost_thresh=cfg.get('tile_cost_thresh',
+                                                            500),
+                                   min_region_size=cfg.get('min_region_size',
+                                                           100),
+                                   tile_edge_weight=cfg.get('tile_edge_weight',
+                                                            2.5),
+                                   secondary_arc_flow_max=cfg.get(
+                                       'secondary_arc_flow_max', 8),
+                                   single_tile_reoptimize=cfg.get(
+                                       'single_tile_reoptimize', False))
+    else:
+        tile = None
+    return tile
+
+
+def set_solver_params(snaphu_cfg):
+    """
+    Set user-defined solver parameter options
+    Parameters
+    ----------
+    snaphu_cfg: dict
+      Dictionary containing SNAPHU options
+    Returns:
+    -------
+    solver: isce3.unwrap.snaphu.SolverParams() or None
+      Object containing solver parameters options or None
+    """
+    # If 'solver_parameters' is in snaphu_cfg, inspect setted
+    # options. If None found, assign default
+    if snaphu_cfg.get('solver_parameters') is not None:
+        cfg = snaphu_cfg['solver_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        solver = snaphu.SolverParams(max_flow_inc=cfg.get('max_flow_inc', 4),
+                                     init_max_flow=cfg.get('initial_max_flow',
+                                                           9999),
+                                     arc_max_flow_const=cfg.get(
+                                         'arc_max_flow_const', 3),
+                                     threshold=cfg.get('threshold', 0.001),
+                                     max_cost=cfg.get('max_cost', 1000.0),
+                                     cost_scale=cfg.get('cost_scale', 100.0),
+                                     n_cycle=cfg.get('n_cycle', 200),
+                                     max_new_node_const=cfg.get(
+                                         'max_new_node_const', 0.0008),
+                                     max_n_flow_cycles=cfg.get(
+                                         'max_n_flow_cycles', None),
+                                     max_cycle_frac=cfg.get('max_cycle_frac',
+                                                            0.00001),
+                                     n_conn_node_min=cfg.get('n_conn_node_min',
+                                                             0),
+                                     n_major_prune=cfg.get('n_major_prune',
+                                                           2000000000),
+                                     prune_cost_thresh=cfg.get(
+                                         'prune_cost_threshold', 2000000000))
+    else:
+        solver = None
+    return solver
+
+
+def set_connected_components_params(snaphu_cfg):
+    """
+    Set user-defined connected components parameter options
+    Parameters
+    ----------
+    snaphu_cfg: dict
+      Dictionary containing SNAPHU options
+    Returns:
+    -------
+    conn: isce3.unwrap.snaphu.ConnCompParams() or None
+      Object containing connected components parameters options or None
+    """
+    # If 'connected_components_parameters' is in snaphu_cfg,
+    # inspect setted options. If None found, assign default
+
+    if snaphu_cfg.get('connected_components_parameters') is not None:
+        cfg = snaphu_cfg['connected_components_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        conn = snaphu.ConnCompParams(
+            min_frac_area=cfg.get('min_frac_area', 0.01),
+            cost_thresh=cfg.get('cost_threshold', 300),
+            max_ncomps=cfg.get('max_ncomps', 32))
+    else:
+        conn = None
+    return conn
+
+
+def set_corr_bias_params(snaphu_cfg):
+    """
+    Set user-defined correlation bias parameter options
+    Parameters
+    ----------
+    snaphu_cfg: dict
+      Dictionary containing SNAPHU options
+    Returns:
+    -------
+    corr: isce3.unwrap.snaphu.CorrBiasModelParams or None
+      Object containing correlation bias model parameters options or None
+    """
+    if snaphu_cfg.get('correlation_bias_parameters') is not None:
+        cfg = snaphu_cfg['correlation_bias_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        corr = snaphu.CorrBiasModelParams(c1=cfg.get('c1', 1.3),
+                                          c2=cfg.get('c2', 0.14),
+                                          min_corr_factor=cfg.get(
+                                              'min_corr_factor', 1.25))
+    else:
+        corr = None
+    return corr
+
+
+def set_phase_std_params(snaphu_cfg):
+    """
+    Set user-defined phase standard deviation
+    model parameter options
+    Parameters
+    ----------
+    snaphu_cfg: dict
+      Dictionary containing SNAPHU options
+     Returns:
+     -------
+    std: isce3.unwrap.snaphu.PhaseStddevModelParams or None
+      Object containing phase std model parameters options or None
+    """
+    if snaphu_cfg.get('phase_std_parameters') is not None:
+        cfg = snaphu_cfg['phase_std_parameters']
+        cfg = clean_none_cfg_keys(cfg)
+        std = snaphu.PhaseStddevModelParams(c1=cfg.get('c1', 0.4),
+                                            c2=cfg.get('c2', 0.35),
+                                            c3=cfg.get('c3', 0.06),
+                                            sigsq_min=cfg.get('sigsq_min', 1))
+    else:
+        std = None
+    return std
+
+
+def compute_nlooks(rg_igram_step, az_igram_step,
+                   rg_res, az_res):
+    '''
+    Compute effective number of looks
+    Parameters
+    ----------
+    rg_igram_step: float
+      Interferogram spacing along slant range
+    az_igram_step: float
+      Interferogram spacing along azimuth
+    rg_res : float
+      Slant range resolution
+    az_res: float
+      Azimuth resolution
+    Returns:
+    -------
+    nlooks: float
+       Effective number of looks
+    '''
+
+    nlooks = (rg_igram_step * az_igram_step) / (rg_res * az_res)
+    return nlooks
+
+
+def unwrap_snaphu(snaphu_cfg, ugram, conn_comp, igram, corr,
+                  nlooks, scratchdir=None, power=None, mask=None, unw_est=None):
+    """
+    Unwrap igram using SNAPHU and user-defined parameters
+    and ancillary rasters
+    Parameters:
+    snaphu_cfg: dict
+      Dictionary containing SNAPHU parameter options
+    ugram: isce3.io.gdal.Raster
+      Unwrapped interferogram raster (gdal.GDT_Float32)
+    conn_comp: isce3.io.gdal.Raster
+      Connected component raster (gdal.GDT_Byte)
+    igram: isce3.io.gdal.Raster
+      Complex wrapped interferogram raster (gdal.GDT_CFloat32)
+    corr: isce3.io.gdal.Raster
+      Normalized interferometric correlation amplitude (gdal.GDT_Float32)
+    nlooks: float
+      Number of effective looks
+    scratchdir: str
+      Path to scratch directory
+    power: isce3.io.gdal.Raster (or None)
+      Slc power raster in the reference SLC geometry
+    mask: isce3.io.gdal.Raster (or None)
+      Binary mask to apply during unwrapping (gdal.GDT_Byte)
+    unw_est: isce3.io.gdal.Raster (or None)
+      Coarse unwrapped phase estimate (gdal.GDT_Float32)
+    """
+    # Get cost mode and corresponding cost parameter object
+    # assign deformation if not found
+    snaphu_cfg = clean_none_cfg_keys(snaphu_cfg)
+    cost_mode = snaphu_cfg.get('cost_mode', 'defo')
+    cost_params = select_cost_options(snaphu_cfg, cost_mode)
+
+    # Get initialization method. If not found, assign MCF
+    init_method = snaphu_cfg.get('initialization_method', 'mcf')
+
+    # Get tiling parameters
+    tile_params = set_tile_params(snaphu_cfg)
+    # Get solver parameters
+    solver_params = set_solver_params(snaphu_cfg)
+    # Get connected components parameters
+    conn_comp_params = set_connected_components_params(snaphu_cfg)
+    # Get correlation bias model parameters
+    corr_model_params = set_corr_bias_params(snaphu_cfg)
+    # Get phase standard deviation model parameters
+    phase_std_params = set_phase_std_params(snaphu_cfg)
+
+    verbose = snaphu_cfg['verbose'] if snaphu_cfg.get(
+        'verbose') is not None else False
+    debug = snaphu_cfg['debug'] if snaphu_cfg.get(
+        'debug') is not None else False
+
+    # All parameters are set. Ready to unwrap
+    snaphu.unwrap(ugram, conn_comp, igram, corr,
+                  nlooks, cost_mode, cost_params, init_method, power,
+                  mask, unw_est, tile_params, solver_params,
+                  conn_comp_params, corr_model_params, phase_std_params,
+                  scratchdir, verbose, debug)
+
+
+def main(cfg):
+    """
+    Main function. cfg: snaphu dictionary
+    """
+
+    # Check that input files are valid
+    check_input_files(cfg)
+
+    # Check output directory
+    if (cfg.get('scratchdir') is not None) and (
+            not os.path.exists(cfg.get('scratchdir'))):
+        os.mkdir(cfg['scratchdir'], exist_ok=True)
+
+    # Loop over interferograms to unwrap
+    for k in range(len(cfg['igram_raster'])):
+        # Put interferogram and correlation in a ISCE3 raster
+        igram_raster = isce3.io.gdal.Raster(cfg['igram_raster'][k])
+        corr_raster = isce3.io.gdal.Raster(cfg['correlation_raster'][k])
+
+        # If ugram_raster and connected components are None use igram basename
+        if cfg['unwrapped_raster'] is None:
+            basename = os.path.basename(cfg['igram_raster'][k])
+            ugram_name = f'{basename}.ugram'
+            conn_comp_name = ugram_name.replace('ugram', 'conn_comp')
+        else:
+            ugram_name = cfg['unwrapped_raster'][k]
+            conn_comp_name = cfg['connected_components_raster'][k]
+
+        # Create unwrapped and connected components rasters.
+        ugram_raster = isce3.io.gdal.Raster(ugram_name, igram_raster.width,
+                                            igram_raster.length, np.float32,
+                                            "ENVI")
+        conn_comp_raster = isce3.io.gdal.Raster(conn_comp_name,
+                                                igram_raster.width,
+                                                igram_raster.length, np.uint32,
+                                                "ENVI")
+
+        # Check allocation of power, mask, unwrapped phase estimate
+        power_raster = get_raster(cfg, 'power', k)
+        mask_raster = get_raster(cfg, 'mask', k)
+        unw_est_raster = get_raster(cfg, 'unwrapped_phase_estimate', k)
+
+        # Proceed to unwrapping
+        unwrap_snaphu(cfg, ugram_raster, conn_comp_raster, igram_raster,
+                      corr_raster, cfg.get('nlooks'), cfg.get('scratchdir'),
+                      power_raster, mask_raster, unw_est_raster)
+
+
+if __name__ == '__main__':
+    # Get command line parser
+    opts = command_line_parser()
+
+    # Load runconfig
+    runconfig = load_runconfig(opts)
+
+    # Main function
+    main(runconfig['snaphu'])