"""InVEST Coastal Vulnerability"""
import time
import os
import math
import logging
import pickle

import numpy
from osgeo import gdal
from osgeo import osr
from osgeo import ogr
import pandas
import rtree
import shapely
import shapely.wkb
import shapely.ops
import shapely.speedups
import shapely.errors
import pygeoprocessing
import taskgraph

from . import utils
from . import validation

LOGGER = logging.getLogger(__name__)

_N_FETCH_RAYS = 16


def execute(args):
    """InVEST Coastal Vulnerability Model.

    For points along a coastline, evaluate the relative exposure of points
    to coastal hazards based on up to eight biophysical hazard indices.
    Also quantify the role of habitats in reducing the hazard. Optionally
    summarize the population density in proximity to each shore point.

    Parameters:
        args['workspace_dir'] (string): (required) a path to the directory that
            will write output and other temporary files during calculation.
        args['results_suffix'] (string): (optional) appended to any output
            filename.
        args['aoi_vector_path'] (string): (required) path to a polygon vector
            that is projected in a coordinate system with units of meters.
            The polygon should intersect the landmass and the shelf contour line.
        args['model_resolution'] (string): (required) distance in meters.
            Points are spaced along the coastline at intervals of this distance.
        args['landmass_vector_path'] (string): (required) path to a polygon
            vector representing landmasses in the region of interest.
        args['wwiii_vector_path'] (string): (required) path to a point vector
            containing wind and wave information across the region of interest.
        args['max_fetch_distance'] (string): (required) maximum distance
            in meters to extend rays from shore points. Points with rays equal
            to this distance will accumulate ocean-driven wave exposure along
            those rays and local-wind-driven wave exposure along the shorter rays.
        args['shelf_contour_vector_path'] (string): (required) path to a
            polyline vector delineating edges of the continental shelf
            or other bathymetry contour.
        args['dem_path'] (string): (required) path to a raster representing the
            elevation on land in the region of interest.
        args['dem_averaging_radius'] (int or float): (required) a value >= 0.
            The radius in meters around each shore point in which to compute
            the average elevation.
        args['habitat_table_path'] (string): (rqeuired) path to a CSV file with
            the following four fields:
                'id': unique string to represent each habitat
                'path': absolute or relative path to a polygon vector
                'rank': integer from 1 to 5 representing the relative
                    protection offered by this habitat
                'protection distance (m)': integer or float used as a
                    search radius around each shore point.
        args['geomorphology_vector_path'] (string): (optional) path to a
            polyline vector that has a field called "RANK" with values from
            1 to 5 in the attribute table.
        args['geomorphology_fill_value'] (int): (optional) a value from 1 to 5
            that will be used as a geomorphology rank for any points not
            proximate to the geomorphology_vector_path.
        args['population_rater_path'] (string): (optional) path a raster with
            values of total population per pixel.
        args['population_radius'] (int or float): (optional) a value >= 0.
            The radius in meters around each shore point in which to compute
            the population density.
        args['slr_vector_path'] (string): (optional) path to point vector
            containing the field `args['slr_field']`.
        args['slr_field'] (string): name of a field in `args['slr_vector_path']`
            containing numeric values.
        args['n_workers'] (int): (optional) The number of worker processes to
            use for processing this model.  If omitted, computation will take
            place in the current process.

    Returns:
        None

    """
    LOGGER.info('Validating arguments')
    invalid_parameters = validate(args)
    if invalid_parameters:
        raise ValueError("Invalid parameters passed: %s" % invalid_parameters)
    _validate_habitat_table_header(args['habitat_table_path'])

    output_dir = os.path.join(args['workspace_dir'])
    intermediate_dir = os.path.join(
        args['workspace_dir'], 'intermediate')
    habitat_dir = os.path.join(
        intermediate_dir, 'habitats')
    shore_dir = os.path.join(
        intermediate_dir, 'shore_points')
    relief_dir = os.path.join(
        intermediate_dir, 'relief')
    geomorph_dir = os.path.join(
        intermediate_dir, 'geomorphology')
    wind_wave_dir = os.path.join(
        intermediate_dir, 'wind_wave')
    surge_dir = os.path.join(
        intermediate_dir, 'surge')
    population_dir = os.path.join(
        intermediate_dir, 'population')
    slr_dir = os.path.join(
        intermediate_dir, 'sealevelrise')

    utils.make_directories(
        [output_dir, intermediate_dir, habitat_dir, shore_dir, relief_dir,
         geomorph_dir, wind_wave_dir, surge_dir, population_dir, slr_dir])
    file_suffix = utils.make_suffix_string(args, 'results_suffix')

    taskgraph_cache_dir = os.path.join(intermediate_dir, '_taskgraph_working_dir')
    try:
        n_workers = int(args['n_workers'])
    except (KeyError, ValueError, TypeError):
        # KeyError when n_workers is not present in args
        # ValueError when n_workers is an empty string.
        # TypeError when n_workers is None.
        n_workers = -1  # Single process mode.
    task_graph = taskgraph.TaskGraph(taskgraph_cache_dir, n_workers)

    model_resolution = float(args['model_resolution'])
    max_fetch_distance = float(args['max_fetch_distance'])

    aoi_vector_info = pygeoprocessing.get_vector_info(
        args['aoi_vector_path'])
    aoi_srs_wkt = aoi_vector_info['projection']
    aoi_bounding_box = aoi_vector_info['bounding_box']
    # add the max_fetch_distance to the bounding box so we can use
    # the clipped landmass in the ray casting routine.
    fetch_buffer = max_fetch_distance + model_resolution
    aoi_bounding_box[0] -= fetch_buffer
    aoi_bounding_box[1] -= fetch_buffer
    aoi_bounding_box[2] += fetch_buffer
    aoi_bounding_box[3] += fetch_buffer

    clipped_landmass_path = os.path.join(
        shore_dir, 'clipped_projected_landmass%s.gpkg' % file_suffix)
    tmp_clipped_path = os.path.join(
        shore_dir, 'tmp_clipped_landmass%s.gpkg' % file_suffix)
    clip_landmass_to_aoi_task = task_graph.add_task(
        func=clip_and_project_vector,
        args=(args['landmass_vector_path'], aoi_bounding_box,
              aoi_srs_wkt, tmp_clipped_path, clipped_landmass_path),
        target_path_list=[clipped_landmass_path, tmp_clipped_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[],
        task_name='clip landmass to aoi')

    target_polygon_pickle_path = os.path.join(
        shore_dir, 'landmass_polygon%s.pickle' % file_suffix)
    target_lines_pickle_path = os.path.join(
        shore_dir, 'landmass_line_index%s.pickle' % file_suffix)
    target_rtree_path = os.path.join(
        shore_dir, 'landmass_line_rtree%s.dat' % file_suffix)
    prep_landmass_geom_task = task_graph.add_task(
        func=prepare_landmass_line_index,
        args=(clipped_landmass_path, target_polygon_pickle_path,
              target_lines_pickle_path, target_rtree_path),
        target_path_list=[target_polygon_pickle_path,
                          target_lines_pickle_path, target_rtree_path],
        ignore_path_list=[target_rtree_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[clip_landmass_to_aoi_task],
        task_name='index landmass geometry')

    shore_point_vector_path = os.path.join(
        shore_dir, 'shore_points%s.gpkg' % file_suffix)
    create_shore_points_task = task_graph.add_task(
        func=interpolate_shore_points,
        args=(args['aoi_vector_path'], target_lines_pickle_path,
              model_resolution, shore_point_vector_path),
        target_path_list=[shore_point_vector_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[
            clip_landmass_to_aoi_task, prep_landmass_geom_task],
        task_name='create shore points')

    exposure_variables_task_list = []
    exposure_variables_path_list = []  # list of 3-tuples like:
    # ('pickle path' (str), bin values (bool), 'R_name' (str))
    # bin values (bool):
        # True when a variable contains values that need binning
            # by percentile to convert to 1-5 ranks.
        # False when variable is already on the 1-5 rank scale.

    target_wwiii_point_vector_path = os.path.join(
        wind_wave_dir, 'wwiii_shore_points%s.gpkg' % file_suffix)
    interpolate_wwiii_task = task_graph.add_task(
        func=interpolate_wwiii_to_shore,
        args=(shore_point_vector_path, args['wwiii_vector_path'],
              target_wwiii_point_vector_path),
        target_path_list=[target_wwiii_point_vector_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[create_shore_points_task],
        task_name='interpolate wwiii to shore points')
    exposure_variables_task_list.append(interpolate_wwiii_task)

    fetch_point_vector_path = os.path.join(
        wind_wave_dir, 'fetch_points%s.gpkg' % file_suffix)
    target_wind_exposure_pickle_path = os.path.join(
        wind_wave_dir, 'wind%s.pickle' % file_suffix)
    exposure_variables_path_list.append(
        (target_wind_exposure_pickle_path, True, 'R_wind'))
    wind_exposure_task = task_graph.add_task(
        func=calculate_wind_exposure,
        args=(target_wwiii_point_vector_path, target_polygon_pickle_path,
              target_rtree_path, target_lines_pickle_path,
              wind_wave_dir, file_suffix, max_fetch_distance,
              fetch_point_vector_path, target_wind_exposure_pickle_path),
        target_path_list=[fetch_point_vector_path,
                          target_wind_exposure_pickle_path],
        ignore_path_list=[target_rtree_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[interpolate_wwiii_task],
        task_name='calculate wind exposure')
    exposure_variables_task_list.append(wind_exposure_task)

    target_wave_exposure_path = os.path.join(
        wind_wave_dir, 'wave%s.pickle' % file_suffix)
    exposure_variables_path_list.append(
        (target_wave_exposure_path, True, 'R_wave'))
    exposure_variables_task_list.append(task_graph.add_task(
        func=calculate_wave_exposure,
        args=(fetch_point_vector_path, max_fetch_distance,
              target_wave_exposure_path),
        target_path_list=[target_wave_exposure_path],
        dependent_task_list=[wind_exposure_task],
        task_name='calculate wave exposure'))

    target_surge_exposure_path = os.path.join(
        surge_dir, 'surge%s.pickle' % file_suffix)
    exposure_variables_path_list.append(
        (target_surge_exposure_path, True, 'R_surge'))
    exposure_variables_task_list.append(task_graph.add_task(
        func=calculate_surge_exposure,
        args=(shore_point_vector_path, args['shelf_contour_vector_path'],
              args['aoi_vector_path'], target_surge_exposure_path),
        target_path_list=[target_surge_exposure_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[create_shore_points_task],
        task_name='calculate surge exposure'))

    relief_point_pickle_path = os.path.join(
        relief_dir, 'relief%s.pickle' % file_suffix)
    exposure_variables_path_list.append(
        (relief_point_pickle_path, True, 'R_relief'))
    exposure_variables_task_list.append(task_graph.add_task(
        func=calculate_relief_exposure,
        args=(shore_point_vector_path, args['dem_path'],
              float(args['dem_averaging_radius']), model_resolution,
              relief_dir, file_suffix, relief_point_pickle_path),
        target_path_list=[relief_point_pickle_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=[create_shore_points_task],
        task_name='calculate relief exposure'))

    # Joining this task instead of passing it to this scheduler function.
    # Tasks added by the scheduler are dependent on shore_points_task
    create_shore_points_task.join()
    hab_tasks_list, hab_targets_list = _schedule_habitat_tasks(
        shore_point_vector_path, args['habitat_table_path'],
        habitat_dir, file_suffix, task_graph)

    target_habitat_protection_path = os.path.join(
        habitat_dir, 'habitat_protection%s.csv' % file_suffix)
    exposure_variables_task_list.append(task_graph.add_task(
        func=calculate_habitat_rank,
        args=(hab_targets_list,
              target_habitat_protection_path),
        target_path_list=[target_habitat_protection_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=hab_tasks_list,
        task_name='calculate habitat protection'))

    if 'geomorphology_vector_path' in args and args['geomorphology_vector_path'] != '':
        projected_geomorphology_vector_path = os.path.join(
            geomorph_dir, 'geomorphology_projected%s.shp' % file_suffix)
        target_srs_wkt = pygeoprocessing.get_vector_info(
            args['aoi_vector_path'])['projection']
        project_geomorph_task = task_graph.add_task(
            func=pygeoprocessing.reproject_vector,
            args=(args['geomorphology_vector_path'], target_srs_wkt,
                  projected_geomorphology_vector_path),
            target_path_list=[projected_geomorphology_vector_path],
            hash_algorithm='md5',
            copy_duplicate_artifact=True,
            task_name='project geomorphology input')

        target_geomorphology_pickle_path = os.path.join(
            geomorph_dir, 'geomorph%s.pickle' % file_suffix)
        exposure_variables_path_list.append(
            (target_geomorphology_pickle_path, False, 'R_geomorph'))
        exposure_variables_task_list.append(task_graph.add_task(
            func=calculate_geomorphology_exposure,
            args=(projected_geomorphology_vector_path,
                  int(args['geomorphology_fill_value']),
                  shore_point_vector_path, model_resolution,
                  target_geomorphology_pickle_path),
            target_path_list=[
                target_geomorphology_pickle_path],
            hash_algorithm='md5',
            copy_duplicate_artifact=True,
            dependent_task_list=[
                create_shore_points_task, project_geomorph_task],
            task_name='calculate geomorphology exposure'))

    if 'slr_vector_path' in args and args['slr_vector_path'] != '':
        target_slr_pickle_path = os.path.join(
            slr_dir, 'slr%s.pickle' % file_suffix)
        exposure_variables_path_list.append(
            (target_slr_pickle_path, True, 'R_slr'))
        exposure_variables_task_list.append(task_graph.add_task(
            func=interpolate_sealevelrise_points,
            args=(shore_point_vector_path, args['slr_vector_path'],
                  args['slr_field'], target_slr_pickle_path),
            target_path_list=[target_slr_pickle_path],
            hash_algorithm='md5',
            copy_duplicate_artifact=True,
            dependent_task_list=[create_shore_points_task],
            task_name='interpolate sea-level rise values'))

    if 'population_raster_path' in args and args['population_raster_path'] != '':
        target_population_pickle_path = os.path.join(
            population_dir, 'population%s.pickle' % file_suffix)
        exposure_variables_path_list.append(
            (target_population_pickle_path, False, None))
        exposure_variables_task_list.append(task_graph.add_task(
            func=aggregate_population_density,
            args=(shore_point_vector_path, args['population_raster_path'],
                  float(args['population_radius']), model_resolution,
                  population_dir, file_suffix, target_population_pickle_path),
            target_path_list=[target_population_pickle_path],
            hash_algorithm='md5',
            copy_duplicate_artifact=True,
            dependent_task_list=[create_shore_points_task],
            task_name='aggregate population raster'))

    target_exposure_vector_path = os.path.join(
        output_dir, 'coastal_exposure%s.gpkg' % file_suffix)
    target_exposure_csv_path = os.path.splitext(
        target_exposure_vector_path)[0] + '.csv'
    target_intermediate_vector_path = os.path.join(
        intermediate_dir, 'intermediate_exposure%s.gpkg' % file_suffix)
    task_graph.add_task(
        func=assemble_results_and_calculate_exposure,
        args=(exposure_variables_path_list, target_habitat_protection_path,
              shore_point_vector_path, target_intermediate_vector_path,
              target_exposure_vector_path, target_exposure_csv_path),
        target_path_list=[
            target_exposure_vector_path, target_exposure_csv_path,
            target_intermediate_vector_path],
        hash_algorithm='md5',
        copy_duplicate_artifact=True,
        dependent_task_list=exposure_variables_task_list,
        task_name='calculate final risk')

    task_graph.close()
    task_graph.join()


def prepare_landmass_line_index(
        landmass_vector_path, target_polygon_pickle_path,
        target_lines_pickle_path, target_rtree_path):
    """Prepare landmass polygon geometry for line operations.

    Converts landmass polygons to lines used in shore point interpolation
    and saves spatial index of lines for fast intersections. Also
    saves unioned geometry polygons for fast point-in-polygon checks.

    Parameters:
        landmass_vector_path (string): path to polygon vector
        target_polygon_pickle_path (string): path to pickle
            storing shapely polygon geometry.
        target_lines_pickle_path (string): path to pickle
            storing list of shapely line geometries
        target_rtree_path (string): path to rtree file indexing
            bounds of line geometries.

    Returns:
        None

    """
    LOGGER.info("preparing landmass geometry")
    # Get shapely geometries from landmass
    landmass_polygon_shapely_list = []
    landmass_polygon_vector = gdal.OpenEx(
        landmass_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    landmass_polygon_layer = landmass_polygon_vector.GetLayer()
    for feature in landmass_polygon_layer:
        geometry = feature.GetGeometryRef()
        shapely_geometry = shapely.wkb.loads(
            geometry.ExportToWkb())
        if shapely_geometry.is_valid:
            landmass_polygon_shapely_list.append(shapely_geometry)
        else:
            LOGGER.warning(
                'landmass FID:%d is excluded due to invalid geometry'
                % feature.GetFID())
        geometry = None
        feature = None
    landmass_polygon_layer = None
    landmass_polygon_vector = None
    landmass_shapely = shapely.ops.cascaded_union(
        landmass_polygon_shapely_list)
    landmass_polygon_shapely_list = None

    # store polygon geom for point-in-poly check later in ray-casting
    with open(target_polygon_pickle_path, 'wb') as polygon_pickle_file:
        pickle.dump(landmass_shapely, polygon_pickle_file)

    # Build index of landmass line geometries on disk
    target_rtree_base = os.path.splitext(target_rtree_path)[0]
    if os.path.exists(target_rtree_path):
        for ext in ['.dat', '.idx']:
            os.remove(target_rtree_base + ext)

    polygon_line_rtree = rtree.index.Index(target_rtree_base)
    shapely_line_index = []
    line_id = 0
    LOGGER.info("indexing geometry of landmass")
    for line in geometry_to_lines(landmass_shapely):
        if (line.bounds[0] == line.bounds[2] and
                line.bounds[1] == line.bounds[3]):
            continue
        polygon_line_rtree.insert(line_id, line.bounds)
        line_id += 1
        shapely_line_index.append(line)

    with open(target_lines_pickle_path, 'wb') as lines_pickle_file:
        pickle.dump(shapely_line_index, lines_pickle_file)
    polygon_line_rtree.close()


def interpolate_shore_points(
        aoi_vector_path, landmass_lines_pickle_path, model_resolution,
        target_vector_path):
    """Create points along a polyline geometry at a specified interval.

    Parameters:
        aoi_vector_path (string): path to polygon vector used to define
            boundaries for adding points.
        landmass_lines_pickle_path (string): path to pickle file
            containing list of shapely linestrings
        model_resolution (float): distance in meters for the point's spacing
        target_vector_path (string): path to .gpkg point vector

    Returns:
        None

    """
    LOGGER.info("creating shore points along edge of landmass")
    # create the spatial reference from the base vector
    aoi_vector_info = pygeoprocessing.get_vector_info(aoi_vector_path)
    aoi_spatial_reference = osr.SpatialReference()
    aoi_spatial_reference.ImportFromWkt(aoi_vector_info['projection'])

    aoi_vector = gdal.OpenEx(aoi_vector_path, gdal.OF_VECTOR)
    aoi_layer = aoi_vector.GetLayer()
    aoi_shapely_list = []
    for feature in aoi_layer:
        geom = feature.GetGeometryRef()
        aoi_shapely_list.append(shapely.wkb.loads(geom.ExportToWkb()))
    aoi_shapely = shapely.ops.cascaded_union(aoi_shapely_list)
    aoi_shapely_prepped = shapely.prepared.prep(aoi_shapely)
    aoi_vector = None
    aoi_layer = None

    with open(landmass_lines_pickle_path, 'rb') as lines_pickle_file:
        shapely_line_list = pickle.load(lines_pickle_file)
    lines_in_aoi_list = []
    for shapely_geom in shapely_line_list:
        # Already clipped landlines to AOI, with max-fetch extent
        # still want to clip lines by actual AOI, for shore point creation
        if aoi_shapely_prepped.intersects(shapely_geom):
            intersected_shapely_geom = aoi_shapely.intersection(shapely_geom)
            if not intersected_shapely_geom.is_empty:
                if intersected_shapely_geom.type == 'MultiLineString':
                    shapely_geom_explode = [
                        shapely.geometry.LineString(x)
                        for x in intersected_shapely_geom]
                    lines_in_aoi_list.extend(shapely_geom_explode)
                if intersected_shapely_geom.type in ['Point', 'MultiPoint']:
                    # a polygon-line intersection could generate a point geom
                    continue
                else:
                    lines_in_aoi_list.append(intersected_shapely_geom)
    unioned_line = shapely.ops.linemerge(lines_in_aoi_list)
    # if none of the lines were disjoint before this linemerge,
    # unioned_line will now be a non-iterable LineString.
    # If some lines were disjoint before the linemerge, unioned_line
    # will now be an iterable MultiLineString.
    try:
        iter(unioned_line)
    except(TypeError):
        unioned_line = [unioned_line]

    point_list = []
    for line in unioned_line:
        n_pts = int(math.ceil(line.length / model_resolution))
        points_along = [
            line.interpolate(float(i) / n_pts, normalized=True)
            for i in range(n_pts)]
        point_list.extend(points_along)

    gpkg_driver = ogr.GetDriverByName("GPKG")
    target_vector = gpkg_driver.CreateDataSource(target_vector_path)
    layer_name = os.path.splitext(
        os.path.basename(target_vector_path))[0].encode('utf-8')
    target_layer = target_vector.CreateLayer(
        layer_name, aoi_spatial_reference, ogr.wkbPoint)

    target_defn = target_layer.GetLayerDefn()
    for point_feature in point_list:
        geometry = ogr.Geometry(ogr.wkbPoint)
        geometry.AddPoint(
            point_feature.coords[0][0], point_feature.coords[0][1])
        feature = ogr.Feature(target_defn)
        feature.SetGeometry(geometry)
        target_layer.CreateFeature(feature)

    target_layer = None
    target_vector = None


def interpolate_wwiii_to_shore(
        base_shore_point_vector_path, wwiii_vector_path,
        target_shore_point_vector_path):
    """Spatial join of Wave Watch 3 data to shore points.

    Finds nearest WW3 points to each shore point and calculates
    a weighted average of values with distance weights.

    Parameters:
        base_shore_point_vector_path (string): path to point vector
        wwiii_vector_path (string): path to point shapefile representing
            the Wave Watch III data.
        target_shore_point_vector_path (string): path to point vector file
            with interpolated wwiii data.

    Returns:
        None

    """
    # Wave Watch III data does not cover the entire planet. Avoid taking values
    # from a WWIII point across the ocean by limiting the search distance.
    # 3 degrees is arbitrary, but seems reasonable and generous given the
    # WWIII points are spaced on a 0.5 degree grid.
    _max_wwiii_distance = 3.0  # degrees

    LOGGER.info("Building spatial index for Wave Watch III points")
    wwiii_rtree = rtree.index.Index()
    wwiii_vector = gdal.OpenEx(wwiii_vector_path, gdal.OF_VECTOR)
    wwiii_layer = wwiii_vector.GetLayer()
    for wwiii_feature in wwiii_layer:
        wwiii_geometry = wwiii_feature.GetGeometryRef()
        wwiii_x = wwiii_geometry.GetX()
        wwiii_y = wwiii_geometry.GetY()
        wwiii_rtree.insert(
            wwiii_feature.GetFID(), (wwiii_x, wwiii_y, wwiii_x, wwiii_y))

    # Copy shore point geometry and create fields for WWIII values
    _copy_point_vector_to_gpkg(
        base_shore_point_vector_path, target_shore_point_vector_path)
    points_vector = gdal.OpenEx(
        target_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    points_layer = points_vector.GetLayer()
    wwiii_defn = wwiii_layer.GetLayerDefn()
    field_names = []
    for field_index in range(wwiii_defn.GetFieldCount()):
        field_defn = wwiii_defn.GetFieldDefn(field_index)
        field_name = field_defn.GetName()
        if field_name in ['I', 'J']:
            continue
        field_names.append(field_name)
        points_layer.CreateField(field_defn)

    wwiii_spatial_reference = osr.SpatialReference()
    wwiii_ref_wkt = pygeoprocessing.get_vector_info(
        wwiii_vector_path)['projection']
    wwiii_spatial_reference.ImportFromWkt(wwiii_ref_wkt)

    if bool(wwiii_spatial_reference.IsProjected()):
        # In case a user decided to project this input on their own,
        # we need a max search distance in meters instead of degrees
        _max_wwiii_distance = 210000  # meters, rough equivalent of 3 degrees

    points_spatial_reference = osr.SpatialReference()
    points_ref_wkt = pygeoprocessing.get_vector_info(
        target_shore_point_vector_path)['projection']
    points_spatial_reference.ImportFromWkt(points_ref_wkt)
    points_to_wwiii_transform = osr.CoordinateTransformation(
        points_spatial_reference, wwiii_spatial_reference)

    LOGGER.info("Interpolating Wave Watch III data to shore points")
    wwiii_field_lookup = {}
    for shore_point_feature in points_layer:
        shore_point_geometry = shore_point_feature.GetGeometryRef()
        # Transform each shore point to match the wwiii SRS
        shore_point_longitude, shore_point_latitude, _ = (
            points_to_wwiii_transform.TransformPoint(
                shore_point_geometry.GetX(), shore_point_geometry.GetY()))
        # From wave watch III points within AOI, get nearest from shore point
        nearest_points = list(wwiii_rtree.nearest(
            (shore_point_longitude, shore_point_latitude,
             shore_point_longitude, shore_point_latitude), 3))[0:3]

        # create placeholders for point geometry and field values
        wwiii_points = numpy.empty((3, 2))
        wwiii_values = numpy.empty((3, len(field_names)))
        for fid_index, fid in enumerate(nearest_points):
            wwiii_feature = wwiii_layer.GetFeature(fid)
            wwiii_geometry = wwiii_feature.GetGeometryRef()
            wwiii_points[fid_index] = numpy.array(
                [wwiii_geometry.GetX(), wwiii_geometry.GetY()])
            try:
                wwiii_values[fid_index] = wwiii_field_lookup[fid]
            except KeyError:
                wwiii_field_lookup[fid] = numpy.array(
                    [float(wwiii_feature.GetField(field_name))
                     for field_name in field_names])
                wwiii_values[fid_index] = wwiii_field_lookup[fid]

        distances = numpy.linalg.norm(
            wwiii_points - numpy.array(
                (shore_point_longitude,
                 shore_point_latitude)), axis=1)

        # make sure points are within a valid data distance
        close_enough = distances < _max_wwiii_distance
        if not any(close_enough):
            raise ValueError(
                'No WaveWatchIII points were found near the area of interest.'
                'Is the area of interest far outside the coverage of %s?'
                % (wwiii_vector_path))

        wwiii_values = numpy.average(
            wwiii_values[close_enough],
            weights=(_max_wwiii_distance - distances[close_enough]),
            axis=0)

        for field_name_index, field_name in enumerate(field_names):
            shore_point_feature.SetField(
                field_name, wwiii_values[field_name_index])
        points_layer.SetFeature(shore_point_feature)
        shore_point_feature = None

    points_layer = None
    points_vector = None


def interpolate_sealevelrise_points(
        base_shore_point_vector_path, slr_points_vector_path,
        slr_fieldname, target_pickle_path):
    """Spatial join of sea-level rise data to shore points.

    Finds nearest sea-level points to each shore point and calculates
    a weighted average of values with inverse-distance weights.

    Parameters:
        base_shore_point_vector_path (string): path to point vector
        slr_points_vector_path (string): path to point vector containing
            the field `slr_fieldname`.
        slr_fieldname (string): name of a field containing numeric values
        target_pickle_path (string): path to pickle file storing dict
            keyed by point fid.

    Returns:
        None

    """
    # In case SLR points map to very far away from the shoreline points
    # due to data preparation error, we don't want to silently join those
    # data to the shoreline points. So, limit the distance to search for SLR
    # points. 500km is arbitrary and generous, but it's only meant to
    # catch errors. A max distance is also useful for inverting distances
    # later on, to use as weights in a weighted average.
    max_slr_distance = 500000  # meters because that unit is mandated for AOI

    base_srs_wkt = pygeoprocessing.get_vector_info(
        slr_points_vector_path)['projection']
    base_spatial_reference = osr.SpatialReference()
    base_spatial_reference.ImportFromWkt(base_srs_wkt)
    target_srs_wkt = pygeoprocessing.get_vector_info(
        base_shore_point_vector_path)['projection']
    target_spatial_reference = osr.SpatialReference()
    target_spatial_reference.ImportFromWkt(target_srs_wkt)
    transform_slr_to_shore = osr.CoordinateTransformation(
        base_spatial_reference, target_spatial_reference)

    slr_rtree = rtree.index.Index()
    slr_vector = gdal.OpenEx(
        slr_points_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    slr_layer = slr_vector.GetLayer()
    for feature in slr_layer:
        geometry = feature.GetGeometryRef()
        geometry.Transform(transform_slr_to_shore)
        slr_x = geometry.GetX()
        slr_y = geometry.GetY()
        slr_rtree.insert(
            feature.GetFID(), (slr_x, slr_y, slr_x, slr_y))
        geometry = None
        feature = None

    result = {}
    base_vector = gdal.OpenEx(
        base_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_layer = base_vector.GetLayer()
    for feature in base_layer:
        shore_fid = feature.GetFID()
        shore_geometry = feature.GetGeometryRef()
        shore_x = shore_geometry.GetX()
        shore_y = shore_geometry.GetY()
        nearest_points = list(slr_rtree.nearest(
            (shore_x, shore_y, shore_x, shore_y), 2))[0:2]

        nearest_slr_values = numpy.empty((len(nearest_points), 1))
        # nearest_slr_values = numpy.empty((2, 1))
        slr_points = numpy.empty((len(nearest_points), 2))
        for fid_index, fid in enumerate(nearest_points):
            slr_feature = slr_layer.GetFeature(fid)
            slr_geometry = slr_feature.GetGeometryRef()
            slr_geometry.Transform(transform_slr_to_shore)
            slr_points[fid_index] = numpy.array(
                [slr_geometry.GetX(), slr_geometry.GetY()])
            try:
                nearest_slr_values[fid_index] = float(
                    slr_feature.GetField(slr_fieldname))
            except KeyError:
                raise KeyError(
                    'fieldname %s not found in %s' % (
                        slr_fieldname, slr_points_vector_path))

        distances = numpy.linalg.norm(
            slr_points - numpy.array(
                (shore_x, shore_y)), axis=1)

        # make sure points are within a valid data distance
        close_enough = distances < max_slr_distance
        if not any(close_enough):
            result[shore_fid] = numpy.nan
            continue

        slr_value = numpy.average(
            nearest_slr_values[close_enough],
            weights=(max_slr_distance - distances[close_enough]),
            axis=0)
        result[shore_fid] = slr_value[0]

    slr_layer = None
    slr_vector = None
    base_layer = None
    base_vector = None

    with open(target_pickle_path, 'wb') as pickle_file:
        pickle.dump(result, pickle_file)


def calculate_wind_exposure(
        base_shore_point_vector_path, landmass_polygon_pickle_path,
        landmass_line_rtree_path, landmass_lines_pickle_path,
        workspace_dir, file_suffix, max_fetch_distance,
        target_shore_point_vector_path, target_wind_exposure_pickle_path):
    """Calculate wind exposure for each shore point.

    Parameters:
        base_shore_point_vector_path (string): path to a point vector
            with WWIII variables in the table.
        landmass_polygon_pickle_path (string): path to pickle
            storing shapely polygon geometry of the landmass.
        landmass_line_rtree_path (string): path to rtree file indexing
            bounds of line geometries.
        landmass_lines_pickle_path (string): path to pickle
            storing list of shapely line geometries. List index must match
            index of `landmass_line_rtree_path`.
        workspace_dir (string): path to a directory for intermediate files
        file_suffix (string): to be appended to output filenames
        max_fetch_distance (float): maximum fetch distance for a ray in
            meters.
        target_shore_point_vector_path (string): path to target point file,
            will be a copy of `base_shore_point_vector_path`'s geometry with
            an 'REI' (relative exposure index) field added.
        target_wind_exposure_pickle_path (string): path to pickle file storing
            dict keyed by shore point fid, with wind exposure values.

    Returns:
        None

    """
    LOGGER.info("Calculating wind exposure")
    temp_fetch_rays_path = os.path.join(
        workspace_dir, 'fetch_rays%s.gpkg' % file_suffix)

    # this should still match the user-defined SRS from the AOI:
    base_ref_wkt = pygeoprocessing.get_vector_info(
        base_shore_point_vector_path)['projection']
    base_spatial_reference = osr.SpatialReference()
    base_spatial_reference.ImportFromWkt(base_ref_wkt)

    gpkg_driver = ogr.GetDriverByName('GPKG')

    # Get prepared polygon geom for point-in-poly checks
    with open(landmass_polygon_pickle_path, 'rb') as polygon_pickle_file:
        landmass_shapely = pickle.load(polygon_pickle_file)
    landmass_shapely_prep = shapely.prepared.prep(landmass_shapely)

    with open(landmass_lines_pickle_path, 'rb') as lines_pickle_file:
        shapely_line_index = pickle.load(lines_pickle_file)

    # load an existing rtree from disk
    polygon_line_rtree = rtree.index.Index(
        os.path.splitext(landmass_line_rtree_path)[0])

    # create fetch rays
    temp_fetch_rays_vector = gpkg_driver.CreateDataSource(
        temp_fetch_rays_path)
    layer_name = os.path.splitext(
        os.path.basename(temp_fetch_rays_path))[0].encode('utf-8')
    temp_fetch_rays_layer = (
        temp_fetch_rays_vector.CreateLayer(
            layer_name, base_spatial_reference, ogr.wkbLineString))
    temp_fetch_rays_defn = temp_fetch_rays_layer.GetLayerDefn()
    temp_fetch_rays_layer.CreateField(ogr.FieldDefn(
        'fetch_dist', ogr.OFTReal))

    base_shore_point_vector = ogr.Open(base_shore_point_vector_path)
    gpkg_driver.CopyDataSource(
        base_shore_point_vector, target_shore_point_vector_path)
    base_shore_point_vector = None

    target_shore_point_vector = gdal.OpenEx(
        target_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    target_shore_point_layer = target_shore_point_vector.GetLayer()
    target_shore_point_layer.CreateField(
        ogr.FieldDefn('REI', ogr.OFTReal))
    for ray_index in range(_N_FETCH_RAYS):
        compass_degree = int(ray_index * 360 / 16.)
        target_shore_point_layer.CreateField(
            ogr.FieldDefn('fdist_%d' % compass_degree, ogr.OFTReal))

    shore_point_logger = _make_logger_callback(
        "Wind exposure %.2f%% complete.", LOGGER)
    # Iterate over every shore point
    LOGGER.info("Casting rays ")
    result_REI = {}
    for shore_point_feature in target_shore_point_layer:
        shore_point_fid = shore_point_feature.GetFID()
        shore_point_logger(
            float(shore_point_fid) /
            target_shore_point_layer.GetFeatureCount())
        rei_value = 0.0
        # Iterate over every ray direction
        for sample_index in range(_N_FETCH_RAYS):
            compass_degree = int(sample_index * 360 / 16.)
            compass_theta = float(sample_index) / _N_FETCH_RAYS * 360
            rei_pct = shore_point_feature.GetField(
                'REI_PCT%d' % int(compass_theta))
            rei_v = shore_point_feature.GetField(
                'REI_V%d' % int(compass_theta))
            cartesian_theta = -(compass_theta - 90)

            # Determine the direction the ray will point
            delta_x = math.cos(cartesian_theta * math.pi / 180)
            delta_y = math.sin(cartesian_theta * math.pi / 180)

            # Start a ray offset from the shore point
            # so that rays start outside of the landmass.
            # Shore points are interpolated onto the coastline,
            # but floating point error results in points being just
            # barely inside/outside the landmass.
            offset = 1  # 1 meter should be plenty
            shore_point_geometry = shore_point_feature.GetGeometryRef()
            point_a_x = (
                shore_point_geometry.GetX() + delta_x * offset)
            point_a_y = (
                shore_point_geometry.GetY() + delta_y * offset)
            point_b_x = point_a_x + delta_x * (
                max_fetch_distance)
            point_b_y = point_a_y + delta_y * (
                max_fetch_distance)
            shore_point_geometry = None

            # build ray geometry so we can intersect it later
            ray_geometry = ogr.Geometry(ogr.wkbLineString)
            ray_geometry.AddPoint(point_a_x, point_a_y)
            ray_geometry.AddPoint(point_b_x, point_b_y)

            # keep a shapely version of the ray so we can do fast intersection
            # with it and the entire landmass
            ray_point_origin_shapely = shapely.geometry.Point(
                point_a_x, point_a_y)

            ray_length = 0.0
            if not landmass_shapely_prep.intersects(
                    ray_point_origin_shapely):
                # the origin is in ocean

                # This algorithm searches for intersections, if one is found
                # the ray updates and a smaller intersection set is determined
                # by experimentation I've found this is significant, but not
                # an order of magnitude, faster than looping through all
                # original possible intersections.  Since this algorithm
                # will be run for a long time, it's worth the additional
                # complexity
                tested_indexes = set()
                while True:
                    intersection = False
                    ray_envelope = ray_geometry.GetEnvelope()
                    for poly_line_index in polygon_line_rtree.intersection(
                            [ray_envelope[i] for i in [0, 2, 1, 3]]):
                        if poly_line_index in tested_indexes:
                            continue
                        tested_indexes.add(poly_line_index)
                        line_shapely = (
                            shapely_line_index[poly_line_index])
                        line_segment = ogr.CreateGeometryFromWkb(line_shapely.wkb)
                        if ray_geometry.Intersects(line_segment):
                            # if the ray intersects the poly line, test if
                            # the intersection is closer than any known
                            # intersection so far
                            intersection_point = ray_geometry.Intersection(
                                line_segment)
                            # offset the dist with model_resolution
                            # update the endpoint of the ray
                            ray_geometry = ogr.Geometry(ogr.wkbLineString)
                            ray_geometry.AddPoint(point_a_x, point_a_y)
                            ray_geometry.AddPoint(
                                intersection_point.GetX(),
                                intersection_point.GetY())
                            intersection = True
                            break
                    if not intersection:
                        break
                # when we get here `min_point` and `ray_length` are the
                # minimum intersection points for the ray and the landmass
                ray_feature = ogr.Feature(temp_fetch_rays_defn)
                ray_length = ray_geometry.Length()
                ray_feature.SetField('fetch_dist', ray_length)
                ray_feature.SetGeometry(ray_geometry)
                temp_fetch_rays_layer.CreateFeature(ray_feature)
            shore_point_feature.SetField(
                'fdist_%d' % compass_degree, ray_length)
            ray_feature = None
            ray_geometry = None
            rei_value += ray_length * rei_pct * rei_v
        shore_point_feature.SetField('REI', rei_value)
        target_shore_point_layer.SetFeature(shore_point_feature)
        result_REI[shore_point_fid] = rei_value

    target_shore_point_layer.SyncToDisk()
    target_shore_point_layer = None
    target_shore_point_vector = None
    temp_fetch_rays_layer.SyncToDisk()
    temp_fetch_rays_layer = None
    temp_fetch_rays_vector = None

    with open(target_wind_exposure_pickle_path, 'wb') as pickle_file:
        pickle.dump(result_REI, pickle_file)


def calculate_wave_exposure(
        base_fetch_point_vector_path, max_fetch_distance,
        target_wave_exposure_pickle_path):
    """Calculate wave exposure values at each shore point.

    Parameters:
        base_fetch_point_vector_path (string): path to a point shapefile that
            contains 16 'WavP_[direction]' fields, 'WavPPCT[direction]'
            fields, 'fdist_[direction]' fields, a single H, and a single T
            field.
        max_fetch_distance (float): max fetch distance before a wind fetch ray
            is terminated.
        target_wave_exposure_pickle_path (string): path to pickle file storing
            dict keyed by shore point fid.

    Returns:
        None

    """
    LOGGER.info("Calculating wave exposure")
    result_Ew = {}
    base_fetch_point_vector = gdal.OpenEx(
        base_fetch_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_fetch_point_layer = base_fetch_point_vector.GetLayer()
    for base_fetch_point_feature in base_fetch_point_layer:
        point_fid = base_fetch_point_feature.GetFID()
        e_local_wedge = (
            0.5 *
            float(base_fetch_point_feature.GetField('H_10PCT'))**2 *
            float(base_fetch_point_feature.GetField('T_10PCT'))) / float(
                _N_FETCH_RAYS)
        e_ocean = 0.0
        e_local = 0.0
        for sample_index in range(_N_FETCH_RAYS):
            compass_degree = int(sample_index * 360 / 16.)
            fdist = base_fetch_point_feature.GetField(
                'fdist_%d' % compass_degree)
            if numpy.isclose(fdist, max_fetch_distance):
                e_ocean += (
                    base_fetch_point_feature.GetField(
                        'WavP_%d' % compass_degree) *
                    base_fetch_point_feature.GetField(
                        'WavPPCT%d' % compass_degree))
            elif fdist > 0.0:
                e_local += e_local_wedge
        result_Ew[point_fid] = max(e_ocean, e_local)

    with open(target_wave_exposure_pickle_path, 'wb') as pickle_file:
        pickle.dump(result_Ew, pickle_file)


def calculate_surge_exposure(
        base_shore_point_vector_path, shelf_contour_path, aoi_vector_path,
        target_surge_pickle_path):
    """Calculate surge potential as distance to nearest point on a contour.

    Parameters:
        base_shore_point_vector_path (string):  path to a point shapefile to
            for relief point analysis.
        shelf_contour_path (string): path to a polyline vector.
        aoi_vector_path (string): path to a polygon vector.
        target_surge_pickle_path (string): path to pickle file storing dict
            keyed by shore point fid.

    Returns:
        None.

    """
    LOGGER.info("Calculating surge potential")
    aoi_vector_info = pygeoprocessing.get_vector_info(
        aoi_vector_path)
    target_srs_wkt = aoi_vector_info['projection']
    target_spatial_reference = osr.SpatialReference()
    target_spatial_reference.ImportFromWkt(target_srs_wkt)
    aoi_bounding_box = aoi_vector_info['bounding_box']

    base_srs_wkt = pygeoprocessing.get_vector_info(
        shelf_contour_path)['projection']
    base_spatial_reference = osr.SpatialReference()
    base_spatial_reference.ImportFromWkt(base_srs_wkt)
    transform = osr.CoordinateTransformation(
        base_spatial_reference, target_spatial_reference)

    base_srs_clipping_box = pygeoprocessing.transform_bounding_box(
        aoi_bounding_box, target_srs_wkt, base_srs_wkt, edge_samples=11)
    base_srs_clipping_geom = ogr.CreateGeometryFromWkt(
        shapely.geometry.box(*base_srs_clipping_box).wkt)

    base_vector = gdal.OpenEx(
        shelf_contour_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_layer = base_vector.GetLayer()

    shapely_geometry_list = []
    for feature in base_layer:
        geometry = feature.GetGeometryRef()
        clipped_geometry = geometry.Intersection(base_srs_clipping_geom)
        if not clipped_geometry.IsEmpty():
            err_code = clipped_geometry.Transform(transform)
            if err_code != 0:
                LOGGER.warning(
                    "Could not transform feature from %s to spatial reference"
                    "system of %s" % (shelf_contour_path, aoi_vector_path))
                continue
            shapely_geom = shapely.wkb.loads(clipped_geometry.ExportToWkb())
            if shapely_geom.type == 'MultiLineString':
                shapely_geom_explode = [
                    shapely.geometry.LineString(x) for x in shapely_geom]
                shapely_geometry_list.extend(shapely_geom_explode)
            else:
                shapely_geometry_list.append(shapely_geom)

    if not shapely_geometry_list:
        raise ValueError(
            "The AOI polygon (%s) and the shelf contour line (%s)"
            "do not intersect" % (aoi_vector_path, shelf_contour_path))

    shelf_shapely_union = shapely.ops.linemerge(
        shapely_geometry_list)

    shore_point_vector = gdal.OpenEx(
        base_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    shore_point_layer = shore_point_vector.GetLayer()
    result = {}
    for point_feature in shore_point_layer:
        point_fid = point_feature.GetFID()
        point_geometry = point_feature.GetGeometryRef()
        point_shapely = shapely.wkb.loads(point_geometry.ExportToWkb())
        result[point_fid] = point_shapely.distance(shelf_shapely_union)

    with open(target_surge_pickle_path, 'wb') as pickle_file:
        pickle.dump(result, pickle_file)


def calculate_relief_exposure(
        base_shore_point_vector_path, base_dem_path, dem_averaging_radius,
        model_resolution, workspace_dir, file_suffix, target_relief_pickle_path):
    """Calculate average of DEM pixels within a radius of shore points.

    Parameters:
        base_shore_point_vector_path (string):  path to a shore point vector.
        base_dem_path (string): path to a DEM raster.
        dem_averaging_radius (float): distance in meters
        model_resolution (float): distance in meters of the shore_point spacing,
            used here as a target pixel size in warp raster.
        workspace_dir (string): path to a directory for intermediate files
        file_suffix (string): to be appended to output filenames
        target_relief_pickle_path (string): path to pickle file storing dict
            keyed by shore point fid.

    Returns:
        None.

    """
    LOGGER.info("Calculating relief exposure")
    base_shore_point_info = pygeoprocessing.get_vector_info(
        base_shore_point_vector_path)
    target_srs_wkt = base_shore_point_info['projection']

    shore_point_bounding_box = base_shore_point_info['bounding_box']
    shore_point_bounding_box[0] -= dem_averaging_radius
    shore_point_bounding_box[1] -= dem_averaging_radius
    shore_point_bounding_box[2] += dem_averaging_radius
    shore_point_bounding_box[3] += dem_averaging_radius

    clipped_projected_dem_path = os.path.join(
        workspace_dir, 'clipped_projected_dem%s.tif' % file_suffix)
    clip_and_project_raster(
        base_dem_path, shore_point_bounding_box, target_srs_wkt,
        model_resolution, workspace_dir, file_suffix,
        clipped_projected_dem_path)

    nodata = pygeoprocessing.get_raster_info(
        clipped_projected_dem_path)['nodata'][0]

    def zero_negative_values(depth_array):
        """Convert negative values to zero."""
        result_array = numpy.empty(depth_array.shape, dtype=numpy.int16)
        if nodata is not None:
            valid_mask = depth_array != nodata
            result_array[:] = nodata
            result_array[valid_mask] = 0
        else:
            result_array[:] = 0
        result_array[depth_array > 0] = depth_array[depth_array > 0]
        return result_array

    positive_dem_path = os.path.join(
        workspace_dir, 'positive_dem%s.tif' % file_suffix)

    pygeoprocessing.raster_calculator(
        [(clipped_projected_dem_path, 1)], zero_negative_values,
        positive_dem_path, gdal.GDT_Int16, nodata)

    _aggregate_raster_values_in_radius(
        base_shore_point_vector_path, positive_dem_path, dem_averaging_radius,
        target_relief_pickle_path, 'mean')


def aggregate_population_density(
        base_shore_point_vector_path, base_population_raster_path,
        search_radius, model_resolution, workspace_dir, file_suffix,
        target_pickle_path):
    """Get population density within a search radius of points.

    Parameters:
        base_shore_point_vector_path (string): path to point vector
        base_population_raster_path (string): path to raster with population
            values
        search_radius (float): radius in meters around each point to search
            for valid population pixels.
        model_resolution (float): distance in meters of the shore_point spacing,
            used here as a target pixel size in warp raster.
        workspace_dir (string): path to directory for intermediate files
        file_suffix (string): appended to any output filename.
        target_pickle_path (string): path to pickle file storing dict
            keyed by shore point fid.

    Returns:
        None

    """
    LOGGER.info("Aggregating population data")
    # SRS here is inherited from the user's AOI
    base_shore_point_info = pygeoprocessing.get_vector_info(
        base_shore_point_vector_path)
    target_srs_wkt = base_shore_point_info['projection']

    # extend the clipping box to accomodate the search radius,
    # plus a little extra to avoid nodata within radius after warping.
    shore_point_bounding_box = base_shore_point_info['bounding_box']
    shore_point_bounding_box[0] -= search_radius
    shore_point_bounding_box[1] -= search_radius
    shore_point_bounding_box[2] += search_radius
    shore_point_bounding_box[3] += search_radius

    clipped_projected_pop_path = os.path.join(
        workspace_dir, 'clipped_projected_pop%s.tif' % file_suffix)
    clip_and_project_raster(
        base_population_raster_path, shore_point_bounding_box, target_srs_wkt,
        model_resolution, workspace_dir, file_suffix, clipped_projected_pop_path)

    _aggregate_raster_values_in_radius(
        base_shore_point_vector_path, clipped_projected_pop_path,
        search_radius, target_pickle_path, 'density')


def _schedule_habitat_tasks(
        base_shore_point_vector_path, habitat_table_path,
        working_dir, file_suffix, task_graph):
    """Add a habitat processing task to the graph, for each habitat.

    Parameters:
        base_shore_point_vector_path (string): path to a shore point vector.
        habitat_table_path (string): path to a CSV file with these fields:
            'id': unique string to represent each habitat
            'path': absolute or relative path to a polygon vector
            'rank': integer from 1 to 5 representing the relative protection
                offered by this habitat
            'protection distance (m)': integer used as a search radius around
                each shore point.
        working_dir (string): path to a directory for intermediate files
        file_suffix (string): to be appended to output filenames
        task_graph (Taskgraph): the graph that was initialized in `execute`

    Returns:
        list of task objects
        list of pickle file path strings

    """
    habitat_dataframe = (pandas.read_csv(
        habitat_table_path, header=0,
        dtype={'id': str, 'path': str, 'rank': numpy.int16,
               'protection distance (m)': numpy.int16})
        .rename(columns={'protection distance (m)': 'distance'}))

    habitat_task_list = []
    habitat_pickles_list = []
    for habitat_row in habitat_dataframe.itertuples():
        base_habitat_path = _sanitize_path(
            habitat_table_path, habitat_row.path)
        target_habitat_pickle_path = os.path.join(
            working_dir, '%s%s.pickle' %
            (habitat_row.id, file_suffix))
        habitat_pickles_list.append(target_habitat_pickle_path)
        habitat_task_list.append(task_graph.add_task(
            func=search_for_habitat,
            args=(base_shore_point_vector_path,
                  habitat_row.distance,
                  habitat_row.rank,
                  habitat_row.id,
                  base_habitat_path,
                  target_habitat_pickle_path),
            target_path_list=[target_habitat_pickle_path],
            task_name='searching for %s' % habitat_row.id))

    return habitat_task_list, habitat_pickles_list


def search_for_habitat(
        base_shore_point_vector_path, search_radius, habitat_rank,
        habitat_id, habitat_vector_path, target_habitat_pickle_path):
    """Search for habitat polygon within a radius of each shore point.

    If habitat is present within the search radius, assign the habitat_rank
    to the shore point FID. If habitat is not present, assign rank of 5.

    Parameters:
        base_shore_point_vector_path (string): path to a shore point vector.
        search_radius (integer): distance around each point to search for
            habitat. units match units from base_shore_point_vector SRS.
        habitat_rank (integer): from 1 to 5 representing the relative
            protection offered by this habitat (5 = no protection).
        habitat_id (string): unique string to represent each habitat.
        habitat_vector_path (string): path to a polygon vector.
        target_habitat_pickle_path (string): path to pickle file storing
            a dict keyed by shore point FID, nested in a dict keyed by
            habitat_id:
                {'habitat_id': {<fid0>: 5, <fid1>.: 5, <fid2>: 5}}

    Returns:
        None

    """
    LOGGER.info("Searching for %s in proximity to shore points" % habitat_id)
    base_shore_point_info = pygeoprocessing.get_vector_info(
        base_shore_point_vector_path)

    shore_point_bounding_box = base_shore_point_info['bounding_box']
    shore_point_bounding_box[0] -= search_radius
    shore_point_bounding_box[1] -= search_radius
    shore_point_bounding_box[2] += search_radius
    shore_point_bounding_box[3] += search_radius

    base_srs_wkt = pygeoprocessing.get_vector_info(
        habitat_vector_path)['projection']
    base_spatial_reference = osr.SpatialReference()
    base_spatial_reference.ImportFromWkt(base_srs_wkt)
    target_srs_wkt = base_shore_point_info['projection']
    target_spatial_reference = osr.SpatialReference()
    target_spatial_reference.ImportFromWkt(target_srs_wkt)
    transform = osr.CoordinateTransformation(
        base_spatial_reference, target_spatial_reference)

    base_srs_clipping_box = pygeoprocessing.transform_bounding_box(
        shore_point_bounding_box, target_srs_wkt, base_srs_wkt, edge_samples=11)
    base_srs_clipping_geom = ogr.CreateGeometryFromWkt(
        shapely.geometry.box(*base_srs_clipping_box).wkt)

    habitat_vector = gdal.OpenEx(
        habitat_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    habitat_layer = habitat_vector.GetLayer()

    # Union all the habitat geometry so that we can easily
    # measure distance from each point to nearest habitat.
    shapely_geometry_list = []
    for feature in habitat_layer:
        geometry = feature.GetGeometryRef()
        clipped_geometry = geometry.Intersection(base_srs_clipping_geom)
        if not clipped_geometry.IsEmpty():
            err_code = clipped_geometry.Transform(transform)
            if err_code != 0:
                LOGGER.warning(
                    "Could not transform feature from %s to"
                    " spatial reference system of AOI"
                    % (habitat_vector_path))
                continue
            shapely_geom = shapely.wkb.loads(clipped_geometry.ExportToWkb())
            if shapely_geom.type == 'MultiPolygon':
                shapely_geom_explode = [
                    shapely.geometry.Polygon(x) for x in shapely_geom]
                shapely_geometry_list.extend(shapely_geom_explode)
            else:
                shapely_geometry_list.append(shapely_geom)
    unioned_habitat_shapely = shapely.ops.cascaded_union(shapely_geometry_list)
    habitat_layer = None
    habitat_vector = None

    base_shore_point_vector = gdal.OpenEx(
        base_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_shore_point_layer = base_shore_point_vector.GetLayer()

    result = {habitat_id: {}}
    for point_feature in base_shore_point_layer:
        point_habitat_rank = 5  # represents no habitat protection
        fid = point_feature.GetFID()

        if unioned_habitat_shapely.is_empty:
            result[habitat_id][fid] = point_habitat_rank
            continue

        point_feature_geometry = point_feature.GetGeometryRef()
        point_shapely = shapely.wkb.loads(
            point_feature_geometry.ExportToWkb())

        # Measure distance to nearest habitat and compare to search radius.
        point_distance_to_feature = point_shapely.distance(
            unioned_habitat_shapely)
        if point_distance_to_feature <= search_radius:
            point_habitat_rank = habitat_rank
        result[habitat_id][fid] = point_habitat_rank

    with open(target_habitat_pickle_path, 'wb') as pickle_file:
        pickle.dump(result, pickle_file)


def calculate_habitat_rank(
        habitat_pickle_list, target_habitat_protection_path):
    """Combine dicts of habitat ranks into a dataframe and calcuate Rhab.

    Parameters:
        habitat_pickle_list (list): list of file paths to pickled dictionaries
            in the form of: {'habitat_id': {<fid0>: 5, <fid1>.: 5, <fid2>: 5}}
        target_habitat_protection_path (string): path to a csv file with a row
            for each shore point, and a header like:
                'fid','kelp','eelgrass','coral','R_hab'

    Returns:
        None

    """
    habitat_id_list = []
    dataframe = pandas.DataFrame()
    for hab_pickle in habitat_pickle_list:
        with open(hab_pickle, 'rb') as file:
            data = pickle.load(file)
        habitat_id = data.keys()[0]  # get habitat name to use in header
        dataframe[habitat_id] = data[habitat_id].values()
        habitat_id_list.append(habitat_id)
    dataframe['fid'] = data[habitat_id].keys()

    def _calc_Rhab(row):
        """Equation 4 from User's Guide."""
        sum_sq_rank = 0.0
        min_rank = 5  # 5 is least protection
        for r in row:
            if r < min_rank:
                min_rank = r
            sum_sq_rank += (5 - r)**2

        if sum_sq_rank > 0:
            r_hab_val = max(
                1.0, 4.8 - 0.5 * (
                    (1.5 * (5 - min_rank))**2 + sum_sq_rank -
                    (5 - min_rank)**2)**0.5)
        else:
            r_hab_val = 5.0
        return r_hab_val

    # Apply _calc_Rhab to each row, excluding the FID value.
    dataframe['R_hab'] = (
        dataframe.drop(columns='fid')
        .apply(axis=1, func=_calc_Rhab))

    # reorder the columns so 'fid' is always first and 'R_hab' last
    (dataframe[['fid'] + habitat_id_list + ['R_hab']]
        .to_csv(target_habitat_protection_path, index=False))


def calculate_geomorphology_exposure(
        geomorphology_vector_path, geomorphology_fill_value,
        base_shore_point_vector_path, model_resolution, target_pickle_path):
    """Join geomorphology ranks to shore points by proximity.

    Buffer each shore point by half the model_resolution and find all
    geomorphology types present around each point. Return the average
    of the geomorphology ranks, unless none are present, then return the

    Parameters:
        geomorphology_vector_path (string): path to polyline vector
            with an integer attribute named `RANK` that contains values
            (1, 2, 3, 4, 5).
        geomorphology_fill_value (int): integer in (1, 2, 3, 4, 5).
        base_shore_point_vector_path (string): path to point vector
        model_resolution (float): distance in meters of the shore_point spacing.
        target_pickle_path (string): path to pickle file storing dict
            keyed by point fid.

    Returns:
        None

    """
    LOGGER.info("Assigning geomorphology rank")
    vector = gdal.OpenEx(geomorphology_vector_path)
    acceptable_ranks = (1, 2, 3, 4, 5)
    geomorph_shapely_lookup = {}
    for rank in acceptable_ranks:
        LOGGER.info("unioning geomorphology rank %d" % rank)
        layer = vector.GetLayer()
        filter_string = ('RANK = %d' % rank)
        layer.SetAttributeFilter(str(filter_string))
        rank_geom_list = []
        for feature in layer:
            geometry = feature.GetGeometryRef()
            shapely_geom = shapely.wkb.loads(geometry.ExportToWkb())
            if not shapely_geom.is_empty:
                if shapely_geom.type == 'MultiLineString':
                    shapely_geom_explode = [
                        shapely.geometry.LineString(x) for x in shapely_geom]
                    rank_geom_list.extend(shapely_geom_explode)
                else:
                    rank_geom_list.append(shapely_geom)
        geomorph_shapely_lookup[rank] = shapely.ops.linemerge(rank_geom_list)
        layer = None
    vector = None

    points_vector = gdal.OpenEx(
        base_shore_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    points_layer = points_vector.GetLayer()
    results = {}
    LOGGER.info("searching for geomorphology in point buffers")
    for point_feature in points_layer:
        fid = point_feature.GetFID()
        point_geom = point_feature.GetGeometryRef()
        poly_geom = point_geom.Buffer(model_resolution / 2.0)
        poly_shapely = shapely.wkb.loads(poly_geom.ExportToWkb())
        geomorph_presence_bool = [
            poly_shapely.intersects(geomorph_shapely_lookup[r]) for r in acceptable_ranks]
        if any(geomorph_presence_bool):
            r_list = [
                r for (r, i) in zip(acceptable_ranks, geomorph_presence_bool) if i]
            mean_geomorph = sum(r_list) / float(len(r_list))
        else:
            mean_geomorph = geomorphology_fill_value
        results[fid] = mean_geomorph

    with open(target_pickle_path, 'wb') as pickle_file:
        pickle.dump(results, pickle_file)


def assemble_results_and_calculate_exposure(
        risk_id_path_list, habitat_protection_path, base_point_vector_path,
        target_intermediate_vector_path, target_output_vector_path,
        target_output_csv_path):
    """Calculate the final exposure score, with and without habitat protection.

    Parameters:
        risk_id_path_list (list): list of 3-tuples like:
            ('relief.pickle', True, 'R_relief')
                1. string: path to pickle with intermediate
                     exposure values for a single variable
                2. bool:
                    if True: variable contains values that need binning
                        by percentile to convert to 1-5 ranks.
                    if False: variable is already on the 1-5 rank scale.
                3. string: if not None, this variable is included in the
                    final exposure score, and 'R_name' becomes a fieldname in
                    target_point_vector_path.
        habitat_protection_path (string): path to csv file with the
           intermediate habitat ranks.
        base_point_vector_path (string): path to shore point vector.
        target_intermediate_vector_path (string): path to point vector
            populated with fields and raw values for all risk variables
            that require binning to ranks.
        target_output_vector_path (string): path to point vector
            populated with fields and rank values for all risk variables,
            final exposure, exposure without habitats, and population density.
        target_output_csv_path (string): path to a csv copy of 
            target_output_vector_path.

    Returns:
        None.

    """
    LOGGER.info("Assembling exposure variables")

    R_hab_name = 'R_hab'

    # To store the final exposure score and variables the score is composed of:
    _copy_point_vector_to_gpkg(
        base_point_vector_path, target_output_vector_path)
    target_output_vector = gdal.OpenEx(
        target_output_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    output_layer = target_output_vector.GetLayer()
    output_layer.CreateField(ogr.FieldDefn(R_hab_name, ogr.OFTReal))

    # To store intermediate values of variables prior to binning into ranks:
    _copy_point_vector_to_gpkg(
        base_point_vector_path, target_intermediate_vector_path)
    target_intermediate_vector = gdal.OpenEx(
        target_intermediate_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    intermediate_layer = target_intermediate_vector.GetLayer()

    # Copy precalculated R_hab to the output vector
    habitat_dataframe = pandas.read_csv(habitat_protection_path)
    for feature_id, rank in zip(
            habitat_dataframe['fid'], habitat_dataframe[R_hab_name]):
        feature = output_layer.GetFeature(int(feature_id))
        feature.SetField(str(R_hab_name), float(rank))
        output_layer.SetFeature(feature)

    # Intermediate variables besides habitat are in pickles.
    # Bin values to ranks if needed, and copy ranks to output vector.
    # and copy pre-binned values to intermediate vector.
    final_values_dict = {}
    intermediate_values_dict = {}
    for pickle_path, to_bins, R_name in risk_id_path_list:
        if R_name:
            fieldname = R_name
        else:
            fieldname = os.path.basename(os.path.splitext(pickle_path)[0])
        exposure_field = ogr.FieldDefn(str(fieldname), ogr.OFTReal)
        exposure_field.SetWidth(24)
        exposure_field.SetPrecision(11)
        output_layer.CreateField(exposure_field)

        if to_bins:
            with open(pickle_path, 'rb') as file:
                raw_values_dict = pickle.load(file)
            intermediate_layer.CreateField(exposure_field)
            invert_values = False
            if R_name == 'R_relief':
                invert_values = True
            values_dict = _bin_values_to_percentiles(
                raw_values_dict, invert_values)
            final_values_dict[fieldname] = values_dict
            intermediate_values_dict[fieldname] = raw_values_dict
        else:
            with open(pickle_path, 'rb') as file:
                final_values_dict[fieldname] = pickle.load(file)

    output_layer.StartTransaction()
    for feature in output_layer:
        fid = feature.GetFID()
        for fieldname in final_values_dict:
            feature.SetField(
                str(fieldname), float(final_values_dict[fieldname][fid]))
        output_layer.SetFeature(feature)
    output_layer.CommitTransaction()
    output_layer = None
    target_output_vector.FlushCache()
    target_output_vector = None

    # Copy pre-binned intermediate values to intermediate output vector
    intermediate_layer.StartTransaction()
    for feature in intermediate_layer:
        fid = feature.GetFID()
        for fieldname in intermediate_values_dict:
            feature.SetField(
                str(fieldname), float(intermediate_values_dict[fieldname][fid]))
        intermediate_layer.SetFeature(feature)
    intermediate_layer.CommitTransaction()
    intermediate_layer = None
    target_intermediate_vector.FlushCache()
    target_intermediate_vector = None

    calculate_final_risk(target_output_vector_path, target_output_csv_path)


def calculate_final_risk(output_vector_path, output_csv_path):
    """Apply geometric mean calculation to variables at each shore point.

    This function modifies the 'output_vector_path' file by adding new fields.

    Parameters:
        output_vector_path (string): path to shore point vector populated
            with numeric fields named with the 'R_' prefix. These fields
            are included in the geometric mean.
        output_csv_path (string): path to csv copy of the final 
            'output_vector_path'.

    Returns:
        None

    """
    LOGGER.info("calculating final risk scores")
    # These fields are added to the output vector and calculated.
    final_risk = 'exposure'
    habitat_role = 'habitat_role'
    risk_no_habitat = 'exposure_no_habitats'

    point_vector = gdal.OpenEx(
        output_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
    point_layer = point_vector.GetLayer()
    point_layer.CreateField(ogr.FieldDefn(final_risk, ogr.OFTReal))
    point_layer.CreateField(ogr.FieldDefn(habitat_role, ogr.OFTReal))
    point_layer.CreateField(ogr.FieldDefn(risk_no_habitat, ogr.OFTReal))
    risk_id_list = [
        field.GetName() for field in point_layer.schema if 'R_' in field.GetName()]

    point_layer.StartTransaction()
    for point_feature in point_layer:
        r_array = numpy.array([
            point_feature.GetField(risk_id)
            for risk_id in risk_id_list], dtype=numpy.float)
        r_tot = _geometric_mean(r_array)
        point_feature.SetField(final_risk, float(r_tot))
        # User's Guide: calculate "coastal_exposure_no_habitats"
        # by replacing R_hab with 5:
        risk_id_list_no_hab = [
            risk_id for risk_id in risk_id_list if risk_id != 'R_hab']
        r_array_nohab = numpy.array([
            point_feature.GetField(risk_id)
            for risk_id in risk_id_list_no_hab] + [5.0], dtype=numpy.float)
        r_tot_no_hab = _geometric_mean(r_array_nohab)
        point_feature.SetField(
            risk_no_habitat, float(r_tot_no_hab))
        point_feature.SetField(
            habitat_role, abs(float(r_tot) - float(r_tot_no_hab)))
        point_layer.SetFeature(point_feature)
    point_layer.CommitTransaction()
    point_vector.FlushCache()

    csv_driver = gdal.GetDriverByName('CSV')
    vector_copy = csv_driver.CreateCopy(
        output_csv_path, point_vector)
    vector_copy = None  # forces flush to disk

    point_layer = None
    point_vector = None


def _geometric_mean(array):
    """Calculate a geometric mean of numpy array of floats.

    Returns:
        float, unless array contains a `nan`, then returns `nan`.

    """
    return numpy.prod(array)**(1./len(array))


def _bin_values_to_percentiles(base_values_dict, invert_values=False):
    """Bin continuous values into percentile categories.

    `nan` should be ignored when calculating percentiles, and features
    with a `nan` value should get a `nan` rank as well.

    Parameters:
        base_values_dict (dictionary): values are numeric
        invert_values (bool): if True, flip sign of values before
            binning.

    Returns:
        Dictionary with same keys as base_values_dict and values
        of integers from 1:5.

    """
    base_values = base_values_dict.values()
    if invert_values:
        base_values = numpy.multiply(base_values, -1.0)

    percentiles = [20, 40, 60, 80, 100]
    rank_array = numpy.searchsorted(
        numpy.nanpercentile(base_values, percentiles),
        base_values).astype(numpy.float) + 1  # floats so we can hold nan
    # nan in base_values need to be re-assigned to nan:
    rank_array[rank_array > len(percentiles)] = numpy.nan

    fids = base_values_dict.keys()
    target_values_dict = {}
    for fid, rank in zip(fids, rank_array):
        target_values_dict[fid] = rank
    return target_values_dict


def clip_and_project_raster(
        base_raster_path, clipping_box, target_srs_wkt, model_resolution,
        working_dir, file_suffix, target_raster_path):
    """Clip a raster to a box in the raster's native SRS, then reproject.

    Parameters:
        base_raster_path (string): path to a gdal raster
        clipping_box (list): sequence of floats that are coordinates in the
            target_srs [minx, miny, maxx, maxy]
        target_srs_wkt (string): well-known-text spatial reference system
        model_resolution (float): value for target pixel size
        working_dir (string): path to directory for intermediate files
        file_suffix (string): appended to any output filename.
        target_raster_path (string): path to clipped and warped raster.

    Returns:
        None

    """
    base_srs_wkt = pygeoprocessing.get_raster_info(
        base_raster_path)['projection']

    # 'base' and 'target' srs are with respect to the base and target raster,
    # so first the clipping box needs to go from 'target' to 'base' srs
    base_srs_clipping_box = pygeoprocessing.transform_bounding_box(
        clipping_box, target_srs_wkt, base_srs_wkt, edge_samples=11)

    clipped_raster_path = os.path.join(
        working_dir,
        os.path.basename(
            os.path.splitext(
                base_raster_path)[0]) + '_clipped%s.tif' % file_suffix)

    base_pixel_size = pygeoprocessing.get_raster_info(
        base_raster_path)['pixel_size']
    target_pixel_size = (model_resolution, model_resolution * -1)

    # Clip in the raster's native srs
    pygeoprocessing.warp_raster(
        base_raster_path, base_pixel_size, clipped_raster_path,
        'bilinear', target_bb=base_srs_clipping_box)

    # Warp to the target SRS
    pygeoprocessing.warp_raster(
        clipped_raster_path, target_pixel_size, target_raster_path,
        'bilinear', target_sr_wkt=target_srs_wkt)


def clip_and_project_vector(
        base_vector_path, clipping_box, target_srs_wkt,
        tmp_vector_path, target_vector_path):
    """Clip a vector to a box in the vector's native SRS, then reproject.

    Parameters:
        base_vector_path (string): path to polygon or multipolygon type vector
        clipping_box (list): sequence of floats that are coordinates in the
            target_srs [minx, miny, maxx, maxy]
        target_srs_wkt (string): well-known-text spatial reference system
        tmp_vector_path (string): path to clipped but unprojected .gpkg vector
        target_vector_path (string): path to clipped and projected .gpkg vector

    Returns:
        None

    """
    for path in [tmp_vector_path, target_vector_path]:
        if os.path.exists(path):
            os.remove(path)

    base_srs_wkt = pygeoprocessing.get_vector_info(
        base_vector_path)['projection']
    base_spatial_reference = osr.SpatialReference()
    base_spatial_reference.ImportFromWkt(base_srs_wkt)

    # 'base' and 'target' srs are with respect to the base and target vector,
    # so first the clipping box needs to go from 'target' to 'base'
    base_srs_clipping_box = pygeoprocessing.transform_bounding_box(
        clipping_box, target_srs_wkt, base_srs_wkt, edge_samples=11)
    base_srs_clipping_box_shapely = shapely.geometry.box(
        *base_srs_clipping_box)

    gpkg_driver = ogr.GetDriverByName("GPKG")
    clipped_vector = gpkg_driver.CreateDataSource(
        tmp_vector_path)
    layer_name = os.path.splitext(
        os.path.basename(tmp_vector_path))[0].encode('utf-8')
    clipped_layer = (
        clipped_vector.CreateLayer(  # intersection could yield MultiPolygon
            layer_name, base_spatial_reference, ogr.wkbMultiPolygon))
    clipped_defn = clipped_layer.GetLayerDefn()

    base_vector = gdal.OpenEx(
        base_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_layer = base_vector.GetLayer()
    for feature in base_layer:
        try:
            feature_shapely = shapely.wkb.loads(
                feature.GetGeometryRef().ExportToWkb())
        except shapely.errors.WKBReadingError:
            feature = None
            LOGGER.warning(
                "Couldn't load a feature in %s", base_vector_path)
            continue
        if base_srs_clipping_box_shapely.intersects(feature_shapely):
            intersection_shapely = base_srs_clipping_box_shapely.intersection(
                feature_shapely)
            clipped_geometry = ogr.CreateGeometryFromWkt(
                intersection_shapely.wkt)
            if clipped_geometry.GetGeometryType() == ogr.wkbPolygon:
                clipped_geometry = ogr.ForceToMultiPolygon(clipped_geometry)
            clipped_feature = ogr.Feature(clipped_defn)
            clipped_feature.SetGeometry(clipped_geometry)
            clipped_layer.CreateFeature(clipped_feature)
            clipped_feature = None

    clipped_layer.SyncToDisk()
    clipped_layer = None
    clipped_vector = None
    base_vector = None
    base_layer = None

    pygeoprocessing.reproject_vector(
        tmp_vector_path, target_srs_wkt, target_vector_path,
        driver_name='GPKG')


def _aggregate_raster_values_in_radius(
        base_point_vector_path, base_raster_path, sample_distance,
        target_pickle_path, aggregation_mode):
    """Aggregate raster values in radius around a point.

    Do the radial search by constructing a rectangular kernel mask
    that approximates a circle.

    E.g. kernel_mask for a raster with pixel_size = (1,1)
    and sample_distance = 3:

    array([[False, False, False,  True, False, False, False],
           [False,  True,  True,  True,  True,  True, False],
           [False,  True,  True,  True,  True,  True, False],
           [ True,  True,  True,  True,  True,  True,  True],
           [False,  True,  True,  True,  True,  True, False],
           [False,  True,  True,  True,  True,  True, False],
           [False, False, False,  True, False, False, False]])

    Parameters:
        base_point_vector_path (string): point vector with projected
            coordinates in units matching sample_disatnce.
        base_raster_path (string): raster file with square pixels and
            projection matching base_point_vector_path.
        sample_distance (float): radius around each point to search
            for valid pixels.
        aggregation_mode (string): 'mean' or 'density'. Calculate the mean
            of valid pixels in search window, or if 'density', divide the mean
            by the area of valid pixels. 'density' assumes units of meters
            for sample_distance and raster geotransform.
        target_pickle_path (string): path to pickle file storing dict
            keyed by point fid.

    Returns:
        None

    """
    if aggregation_mode not in ['mean', 'density']:
        raise ValueError('aggregation mode must be either "mean" or "density"')

    raster = gdal.OpenEx(base_raster_path, gdal.OF_RASTER | gdal.GA_Update)
    band = raster.GetRasterBand(1)
    n_rows = band.YSize
    n_cols = band.XSize
    geotransform = raster.GetGeoTransform()
    nodata = band.GetNoDataValue()

    # we can assume square pixels at this point because
    # we already warped input raster and defined square pixels
    pixel_dist = int(abs(
        sample_distance / (geotransform[1])))

    # kernel dimensions will be 2 * pixel_dist + 1 so that
    # point feature is always inside the center pixel of the kernel.
    # create rectangular kernel and a mask so it looks circular.
    X, Y = numpy.ogrid[:(1 + pixel_dist), :(1 + pixel_dist)]
    lr_quadrant = numpy.hypot(X, Y)
    ll_quadrant = numpy.flip(lr_quadrant[:, 1:], axis=1)
    bottom_half = numpy.concatenate((ll_quadrant, lr_quadrant), axis=1)
    top_half = numpy.flip(bottom_half[1:, :], axis=0)
    kernel_index_distances = numpy.concatenate((top_half, bottom_half), axis=0)

    kernel_mask = numpy.where(
        kernel_index_distances > pixel_dist, False, True)

    result = {}
    vector = gdal.OpenEx(
        base_point_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    layer = vector.GetLayer()
    for point_feature in layer:
        warn_msg = 0
        win_xsize = kernel_mask.shape[0]
        win_ysize = kernel_mask.shape[1]
        temp_kernel_mask = kernel_mask[:]

        point_fid = point_feature.GetFID()
        point_geometry = point_feature.GetGeometryRef()
        point_x = point_geometry.GetX()
        point_y = point_geometry.GetY()
        point_geometry = None

        # kernel origin (upper-left) in pixel-space
        pixel_x = int(
            (point_x - geotransform[0]) /
            geotransform[1]) - pixel_dist
        pixel_y = int(
            (point_y - geotransform[3]) /
            geotransform[5]) - pixel_dist

        # if kernel origin is outside bounds of the raster
        if pixel_x < 0:
            warn_msg = 1
            win_xsize += pixel_x
            pixel_x = 0
            # trim mask columns by indexing backward from the right
            temp_kernel_mask = temp_kernel_mask[:, (-1 * win_xsize):]
        if pixel_y < 0:
            warn_msg = 1
            win_ysize += pixel_y
            pixel_y = 0
            # trim mask rows by indexing backward from the bottom
            temp_kernel_mask = temp_kernel_mask[(-1 * win_ysize):, :]

        # if kernel extent is outside bounds of the raster
        if pixel_x + win_xsize > n_cols:
            warn_msg = 1
            win_xsize -= pixel_x + win_xsize - n_cols
            # trim mask columns by indexing forward from left
            temp_kernel_mask = temp_kernel_mask[:, :win_xsize]
        if pixel_y + win_ysize > n_rows:
            warn_msg = 1
            win_ysize -= pixel_y + win_ysize - n_rows
            # trim mask rows by indexing forward from top
            temp_kernel_mask = temp_kernel_mask[:win_ysize, :]
        if warn_msg:
            LOGGER.warning(
                'search radius around point (%.6f, %.6f) '
                'extends beyond bounds of raster %s'
                % (point_x, point_y, base_raster_path))
        if win_xsize <= 0 or win_ysize <= 0:
            pixel_value = 0
        else:
            try:
                array = band.ReadAsArray(
                    xoff=pixel_x, yoff=pixel_y, win_xsize=win_xsize,
                    win_ysize=win_ysize)
                if nodata is not None:
                    mask = (array != nodata) & temp_kernel_mask
                else:
                    mask = kernel_mask
                if numpy.count_nonzero(mask) > 0:
                    pixel_value = numpy.mean(array[mask])
                else:
                    pixel_value = numpy.nan
            except Exception:
                LOGGER.error(
                    'band size %d %d', band.XSize,
                    band.YSize)
                raise

        # calculate pixel area in sq km
        if aggregation_mode == 'density' and ~numpy.isnan(pixel_value):
            pixel_area_km = abs(
                geotransform[1] * geotransform[5]) / 1e6  # converts m^2 to km^2
            pixel_value /= pixel_area_km

        result[point_fid] = pixel_value
    with open(target_pickle_path, 'wb') as pickle_file:
        pickle.dump(result, pickle_file)


def geometry_to_lines(geometry):
    """Convert a geometry object to a list of lines."""
    if geometry.type == 'Polygon':
        return polygon_to_lines(geometry)
    elif geometry.type == 'MultiPolygon':
        line_list = []
        for geom in geometry.geoms:
            line_list.extend(geometry_to_lines(geom))
        return line_list
    else:
        return []


def polygon_to_lines(geometry):
    """Return a list of shapely lines given higher order shapely geometry."""
    line_list = []
    last_point = geometry.exterior.coords[0]
    for point in geometry.exterior.coords[1::]:
        if point == last_point:
            continue
        line_list.append(shapely.geometry.LineString([last_point, point]))
        last_point = point
    line_list.append(shapely.geometry.LineString([
        last_point, geometry.exterior.coords[0]]))
    for interior in geometry.interiors:
        last_point = interior.coords[0]
        for point in interior.coords[1::]:
            if point == last_point:
                continue
            line_list.append(shapely.geometry.LineString([last_point, point]))
            last_point = point
        line_list.append(shapely.geometry.LineString([
            last_point, interior.coords[0]]))
    return line_list


def _copy_point_vector_to_gpkg(base_vector_path, target_gpkg_path):
    """Copy only the geometry of one GDAL point vector to a geopackage.

    Parameters:
        base_vector_path (string): path to a GDAL-supported vector of
            type wkbPoint.
        target_gpkg_path (string): file path with a '.gpkg' extension.

    Returns:
        None

    """
    gpkg_driver = ogr.GetDriverByName("GPKG")
    aoi_srs = osr.SpatialReference()
    aoi_srs.ImportFromWkt(
        pygeoprocessing.get_vector_info(base_vector_path)['projection'])

    target_vector = gpkg_driver.CreateDataSource(
        target_gpkg_path)
    layer_name = os.path.splitext(
        os.path.basename(target_gpkg_path))[0].encode('utf-8')
    target_layer = target_vector.CreateLayer(
        layer_name, aoi_srs, ogr.wkbPoint)

    target_layer_defn = target_layer.GetLayerDefn()

    base_vector = gdal.OpenEx(base_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
    base_layer = base_vector.GetLayer()
    for feature in base_layer:
        target_geometry = ogr.Geometry(ogr.wkbPoint)
        base_geometry = feature.GetGeometryRef()
        target_geometry.AddPoint(base_geometry.GetX(), base_geometry.GetY())
        target_feature = ogr.Feature(target_layer_defn)
        target_feature.SetGeometry(target_geometry)
        target_layer.CreateFeature(target_feature)
        feature = None
        target_feature = None

    target_layer_defn = None
    target_layer = None
    target_vector = None
    base_layer = None
    base_vector = None


def _sanitize_path(base_path, raw_path):
    """Return `raw_path` if absolute, or make absolute relative to `base_path`."""
    if os.path.isabs(raw_path):
        return raw_path
    return os.path.join(os.path.dirname(base_path), raw_path)


def _make_logger_callback(message, logger):
    """Build a timed logger callback that prints `message` replaced.

    Parameters:
        message (string): a string that expects a %f replacement variable for
            `proportion_complete`.

    Returns:
        Function with signature:
            logger_callback(proportion_complete, psz_message, p_progress_arg)

    """
    def logger_callback(proportion_complete):
        """Argument names come from the GDAL API for callbacks."""
        try:
            current_time = time.time()
            if ((current_time - logger_callback.last_time) > 5.0 or
                    (proportion_complete == 1.0 and
                     logger_callback.total_time >= 5.0)):
                logger.info(message, proportion_complete * 100)
                logger_callback.last_time = current_time
                logger_callback.total_time += current_time
        except AttributeError:
            logger_callback.last_time = time.time()
            logger_callback.total_time = 0.0

    return logger_callback


def _validate_habitat_table_header(habitat_table_path):
    """Validate header row of csv for expected names."""
    dataframe = pandas.read_csv(habitat_table_path)
    header = ('id', 'path', 'rank', 'protection distance (m)')
    missing_names = list(set(header).difference(list(dataframe)))
    if missing_names:
        raise ValueError(
            '%s does not have the following required column names: %s'
            % (habitat_table_path, str(missing_names)))


@validation.invest_validator
def validate(args, limit_to=None):
    """Validate args to ensure they conform to ``execute``'s contract.

    Parameters:
        args (dict): dictionary of key(str)/value pairs where keys and
            values are specified in ``execute`` docstring.
        limit_to (str): (optional) if not None indicates that validation
            should only occur on the ``args[limit_to]`` value. The intent that
            individual key validation could be significantly less expensive
            than validating the entire ``args`` dictionary.

    Returns:
        list of ([invalid key_a, invalid_key_b, ...], 'warning/error message')
            tuples. Where an entry indicates that the invalid keys caused
            the error message in the second part of the tuple. This should
            be an empty list if validation succeeds.

    """
    missing_key_list = []
    no_value_list = []
    validation_error_list = []

    required_keys = [
        'workspace_dir',
        'aoi_vector_path',
        'model_resolution',
        'landmass_vector_path',
        'wwiii_vector_path',
        'max_fetch_distance',
        'shelf_contour_vector_path',
        'dem_path',
        'dem_averaging_radius',
        'habitat_table_path',
        ]

    spatial_files = [
        ('aoi_vector_path', gdal.OF_VECTOR, 'vector'),
        ('landmass_vector_path', gdal.OF_VECTOR, 'vector'),
        ('wwiii_vector_path', gdal.OF_VECTOR, 'vector'),
        ('shelf_contour_vector_path', gdal.OF_VECTOR, 'vector'),
        ('dem_path', gdal.OF_RASTER, 'raster',),
        ]

    polyline_files = [
        'shelf_contour_vector_path'
    ]

    numeric_keys = [
        'max_fetch_distance',
        'dem_averaging_radius',
    ]

    # Some fields are required only when others have values
    if limit_to in [
            None, 'geomorphology_fill_value', 'geomorphology_vector_path']:
        if ('geomorphology_vector_path' in args and
                args['geomorphology_vector_path'] != ''):
            required_keys.append('geomorphology_fill_value')
            spatial_files.append(
                ('geomorphology_vector_path', gdal.OF_VECTOR, 'vector'))
            polyline_files.append('geomorphology_vector_path')
            numeric_keys.append('geomorphology_fill_value')

    if limit_to in [
            None, 'population_raster_path', 'population_radius']:
        if ('population_raster_path' in args and
                args['population_raster_path'] != ''):
            required_keys.append('population_radius')
            spatial_files.append(
                ('population_raster_path', gdal.OF_RASTER, 'raster'))
            numeric_keys.append('population_radius')

    if limit_to in [
            None, 'slr_vector_path', 'slr_field']:
        if ('slr_vector_path' in args and
                args['slr_vector_path'] != ''):
            required_keys.append('slr_field')
            spatial_files.append(
                ('slr_vector_path', gdal.OF_VECTOR, 'vector'))

    for key in required_keys:
        if limit_to in (None, key):
            if key not in args:
                missing_key_list.append(key)
            elif args[key] in ('', None):
                no_value_list.append(key)

    if missing_key_list:
        raise KeyError(*missing_key_list)

    if no_value_list:
        validation_error_list.append(
            (no_value_list, 'parameter has no value'))

    with utils.capture_gdal_logging():
        for key, filetype, filetype_string in spatial_files:
            if limit_to not in (key, None):
                continue

            spatial_file = gdal.OpenEx(args[key], filetype)
            if spatial_file is None:
                validation_error_list.append(
                    ([key], 'Must be a %s' % filetype_string))
                continue

            if key in polyline_files:
                layer = spatial_file.GetLayer()
                # linestring codes: https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogr_core.h
                if layer.GetGeomType() not in [
                        2, 5, 2002, 2005, 3002, 3005, -2147483646, -2147483643]:
                    validation_error_list.append(
                        ([key], 'Must be a polyline vector'))
            if key in 'slr_vector_path':
                layer = spatial_file.GetLayer()
                fields = [field.GetName() for field in layer.schema]
                if args['slr_field'] not in fields:
                    validation_error_list.append(
                        (['slr_field'], 'Must be a fieldname in %s' % args[key]))

    # Verify that the AOI projection is in meters.
    # We don't care about the other spatial inputs, as they'll all be
    # reprojected to the AOI's projection.
    if limit_to in ('aoi_vector_path', None):
        # only do this check if we can open the vector.
        with utils.capture_gdal_logging():
            do_spatial_check = False
            if gdal.OpenEx(args['aoi_vector_path'], gdal.OF_VECTOR) is not None:
                do_spatial_check = True
        if do_spatial_check:
            aoi_srs = osr.SpatialReference()
            aoi_srs.ImportFromWkt(
                pygeoprocessing.get_vector_info(
                    args['aoi_vector_path'])['projection'])
            if (abs(aoi_srs.GetLinearUnits() - 1.0) > 0.5e-7 or
                    not bool(aoi_srs.IsProjected())):
                validation_error_list.append(
                    (['aoi_vector_path'], 'Must be projected in meters'))

    if (limit_to in [None, 'habitat_table_path'] and
            'habitat_table_path' in args):
        if (args['habitat_table_path'] not in [None, ''] and
                not os.path.exists(args['habitat_table_path'])):
            validation_error_list.append(
                (['habitat_table_path'],
                 'not found on disk'))

    for key in numeric_keys:
        if limit_to not in (key, None):
            continue
        try:
            float(args[key])
        except (ValueError, TypeError):
            validation_error_list.append(
                ([key], "Must be a number"))

    return validation_error_list
