"""Module for Regression Testing the InVEST Coastal Vulnerability module."""
import unittest
import tempfile
import shutil
import os
import pickle

from osgeo import gdal, osr, ogr
import numpy.testing
import pandas.testing
from pygeoprocessing.testing import sampledata
from shapely.geometry import Point, Polygon
import taskgraph

from natcap.invest import coastal_vulnerability

REGRESSION_DATA = os.path.join(
    os.path.dirname(__file__), '..', 'data', 'invest-test-data',
    'coastal_vulnerability')
INPUT_DATA = os.path.join(
    os.path.dirname(__file__), '..', 'data', 'invest-test-data',
    'coastal_vulnerability', 'input')


class CoastalVulnerabilityTests(unittest.TestCase):
    """Tests for the Coastal Vulnerability Model."""

    def setUp(self):
        """Overriding setUp function to create temp workspace directory."""
        # this lets us delete the workspace after its done no matter the
        # the rest result
        self.workspace_dir = tempfile.mkdtemp()

    def tearDown(self):
        """Overriding tearDown function to remove temporary directory."""
        shutil.rmtree(self.workspace_dir)

    @staticmethod
    def generate_base_args(workspace_dir):
        """Generate an args dict with default required args."""
        args = {'workspace_dir': workspace_dir,
                'n_workers': -1,
                'wwiii_vector_path': os.path.join(
                    INPUT_DATA, 'WaveWatchIII_subset.shp'),
                'landmass_vector_path': os.path.join(
                    INPUT_DATA, 'land_polygon_simple_utm.shp'),
                'aoi_vector_path': os.path.join(
                    INPUT_DATA, 'AOI_BarkClay.shp'),
                'model_resolution': 25000,
                'max_fetch_distance': 12000,
                'dem_path': os.path.join(
                    INPUT_DATA, 'dem_wgs84.tif'),
                'dem_averaging_radius': 33000.0,
                'habitat_table_path': os.path.join(
                    INPUT_DATA, "natural_habitats_wcvi.csv"),
                'shelf_contour_vector_path': os.path.join(
                    INPUT_DATA, 'continental_shelf_contour.gpkg')
                }
        return args

    def test_wind_and_wave_exposure(self):
        """CV: regression test for wind and wave exposure values.

        The wave calculation depends on an intermediate product from
        the wind calculation, so I'm testing them both in this scope."""

        workspace_dir = self.workspace_dir
        # these points have the WWIII values interpolated onto them.
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, "wwiii_shore_points_5000m.gpkg")
        landmass_polygon_pickle_path = os.path.join(
            INPUT_DATA, "landmass_polygon_5000m.pickle")
        landmass_line_rtree_path = os.path.join(
            INPUT_DATA, "landmass_line_rtree_5000m.dat")
        landmass_lines_pickle_path = os.path.join(
            INPUT_DATA, "landmass_line_index_5000m.pickle")
        file_suffix = ''
        max_fetch_distance = 12000
        target_fetch_point_vector_path = os.path.join(
            workspace_dir, 'fetch_points.gpkg')
        target_wind_exposure_pickle_path = os.path.join(
            workspace_dir, 'wind.pickle')

        coastal_vulnerability.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_fetch_point_vector_path,
            target_wind_exposure_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_wind.pickle')
        assert_pickled_arrays_almost_equal(
            target_wind_exposure_pickle_path, expected_raw_values_path)

        target_wave_exposure_pickle_path = os.path.join(
            workspace_dir, 'wave.pickle')
        coastal_vulnerability.calculate_wave_exposure(
            target_fetch_point_vector_path, max_fetch_distance,
            target_wave_exposure_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_wave.pickle')
        assert_pickled_arrays_almost_equal(
            target_wave_exposure_pickle_path, expected_raw_values_path)

    def test_habitat_rank(self):
        """CV: regression test for habitat ranks."""

        workspace_dir = self.workspace_dir
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, "wwiii_shore_points_5000m.gpkg")
        habitat_table_path = os.path.join(
            INPUT_DATA, "natural_habitats_wcvi.csv")
        target_habitat_protection_path = os.path.join(
            workspace_dir, 'habitat_protection.csv')
        file_suffix = ''

        task_graph = taskgraph.TaskGraph(
            os.path.join(workspace_dir, 'taskgraph_dir'), -1)

        task_list, pickle_list = coastal_vulnerability._schedule_habitat_tasks(
            base_shore_point_vector_path, habitat_table_path,
            workspace_dir, file_suffix, task_graph)

        coastal_vulnerability.calculate_habitat_rank(
            pickle_list, target_habitat_protection_path)

        expected_habitat_path = os.path.join(
            REGRESSION_DATA, 'expected_habitat_protection.csv')
        actual_values_df = pandas.read_csv(target_habitat_protection_path)
        expected_values_df = pandas.read_csv(expected_habitat_path)
        pandas.testing.assert_frame_equal(actual_values_df, expected_values_df)

    def test_geomorphology_rank(self):
        """CV: regression test for geomorphology values."""

        workspace_dir = self.workspace_dir
        geomorphology_vector_path = os.path.join(
            INPUT_DATA, "geomorphology_few_ranks.shp")
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, "wwiii_shore_points_5000m.gpkg")
        target_pickle_path = os.path.join(
            workspace_dir, 'geomorphology.pickle')
        model_resolution = 5000

        coastal_vulnerability.calculate_geomorphology_exposure(
            geomorphology_vector_path, 3,
            base_shore_point_vector_path,
            model_resolution, target_pickle_path)

        expected_values_pickle_path = os.path.join(
            REGRESSION_DATA, 'expected_geomorphology.pickle')
        assert_pickled_arrays_almost_equal(
            target_pickle_path, expected_values_pickle_path)

    def test_surge_exposure_rank(self):
        """CV: regression test for surge exposure values."""

        workspace_dir = self.workspace_dir
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, "wwiii_shore_points_5000m.gpkg")
        shelf_contour_path = os.path.join(
            INPUT_DATA, 'continental_shelf_contour.gpkg')
        aoi_vector_path = os.path.join(
            INPUT_DATA, 'AOI_BarkClay.shp')
        target_surge_pickle_path = os.path.join(
            workspace_dir, 'surge.pickle')

        coastal_vulnerability.calculate_surge_exposure(
            base_shore_point_vector_path, shelf_contour_path,
            aoi_vector_path, target_surge_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_surge.pickle')

        assert_pickled_arrays_almost_equal(
            target_surge_pickle_path, expected_raw_values_path)

    def test_relief_values(self):
        """CV: regression test for aggregated relief values."""

        workspace_dir = self.workspace_dir
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, 'wwiii_shore_points_5000m.gpkg')
        dem_path = os.path.join(
            INPUT_DATA, 'dem_wgs84.tif')
        target_relief_pickle_path = os.path.join(
            workspace_dir, 'relief.pickle')
        dem_averaging_radius = 20000.0
        model_resolution = 5000.0
        file_suffix = ''

        coastal_vulnerability.calculate_relief_exposure(
            base_shore_point_vector_path, dem_path, dem_averaging_radius,
            model_resolution, workspace_dir, file_suffix,
            target_relief_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_relief.pickle')

        # I found minor gdal version differences produced slightly different
        # pixel values from gdal.Warp. So asserting exact values resulting
        # from calculate_relief_exposure might fail depending on gdal version.
        # Rather than pin this test to a specific gdal version, I'm asserting
        # equal rank values after binning the raw values into 1-5 ranks.
        with open(target_relief_pickle_path, 'rb') as file:
            actual_values_dict = pickle.load(file)
        actual_rank_dict = coastal_vulnerability._bin_values_to_percentiles(
            actual_values_dict, invert_values=True)
        with open(expected_raw_values_path, 'rb') as file:
            expected_values_dict = pickle.load(file)
        expected_rank_dict = coastal_vulnerability._bin_values_to_percentiles(
            expected_values_dict, invert_values=True)

        numpy.testing.assert_array_equal(
            actual_rank_dict.values(), expected_rank_dict.values())

    def test_population_values(self):
        """CV: regression test for aggregated population density."""

        workspace_dir = self.workspace_dir
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, 'wwiii_shore_points_5000m.gpkg')
        population_path = os.path.join(
            INPUT_DATA, 'population.tif')
        target_relief_pickle_path = os.path.join(
            workspace_dir, 'population.pickle')
        search_radius = 20000.0
        model_resolution = 5000.0
        file_suffix = ''

        coastal_vulnerability.aggregate_population_density(
            base_shore_point_vector_path, population_path, search_radius,
            model_resolution, workspace_dir, file_suffix,
            target_relief_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_population.pickle')

        assert_pickled_arrays_almost_equal(
            target_relief_pickle_path, expected_raw_values_path)

    def test_interpolate_slr(self):
        """CV: regression test for sea-level rise values.

        This tests an edge case where there is only one point in the
        SLR dataset, and it requires a coordinate transformation.
        """

        workspace_dir = self.workspace_dir
        base_shore_point_vector_path = os.path.join(
            INPUT_DATA, 'wwiii_shore_points_5000m.gpkg')

        # Make an SLR point vector
        slr_fieldname = 'Trend'
        slr_point_vector_path = os.path.join(workspace_dir, 'simple_points.shp')
        out_driver = ogr.GetDriverByName('ESRI Shapefile')
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)
        shapely_feature = Point(-125.65, 49.0)
        out_vector = out_driver.CreateDataSource(slr_point_vector_path)
        layer_name = os.path.basename(os.path.splitext(slr_point_vector_path)[0])
        out_layer = out_vector.CreateLayer(layer_name, srs=srs)
        field_defn = ogr.FieldDefn(slr_fieldname, ogr.OFTReal)
        out_layer.CreateField(field_defn)
        layer_defn = out_layer.GetLayerDefn()
        new_feature = ogr.Feature(layer_defn)
        new_geometry = ogr.CreateGeometryFromWkb(shapely_feature.wkb)
        new_feature.SetGeometry(new_geometry)
        new_feature.SetField(slr_fieldname, 1.3)
        out_layer.CreateFeature(new_feature)
        out_layer = None
        out_vector = None

        target_pickle_path = os.path.join(
            workspace_dir, 'slr.pickle')
        coastal_vulnerability.interpolate_sealevelrise_points(
            base_shore_point_vector_path, slr_point_vector_path,
            slr_fieldname, target_pickle_path)

        expected_raw_values_path = os.path.join(
            REGRESSION_DATA, 'expected_slr.pickle')

        assert_pickled_arrays_almost_equal(
            target_pickle_path, expected_raw_values_path)

    def test_interpolate_slr_beyond_maxdistance(self):
        """CV: test sea-level rise returns nan beyond max search distance."""

        workspace_dir = self.workspace_dir

        # Make a shore point
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N
        shore_point_path = os.path.join(workspace_dir, 'shore_point.shp')
        sampledata.create_vector_on_disk(
            [Point(0., 0.)], srs.ExportToWkt(), filename=shore_point_path)

        # Make an SLR point very far away from the shore point
        slr_fieldname = 'Trend'
        slr_point_vector_path = os.path.join(workspace_dir, 'slr_point.shp')
        out_driver = ogr.GetDriverByName('ESRI Shapefile')
        shapely_feature = Point(1e6, 1e6)  # very far from Point(0, 0)
        out_vector = out_driver.CreateDataSource(slr_point_vector_path)
        layer_name = os.path.basename(os.path.splitext(slr_point_vector_path)[0])
        out_layer = out_vector.CreateLayer(layer_name, srs=srs)
        field_defn = ogr.FieldDefn(slr_fieldname, ogr.OFTReal)
        out_layer.CreateField(field_defn)
        layer_defn = out_layer.GetLayerDefn()
        new_feature = ogr.Feature(layer_defn)
        new_geometry = ogr.CreateGeometryFromWkb(shapely_feature.wkb)
        new_feature.SetGeometry(new_geometry)
        new_feature.SetField(slr_fieldname, 1.3)
        out_layer.CreateFeature(new_feature)
        out_layer = None
        out_vector = None

        target_pickle_path = os.path.join(
            workspace_dir, 'slr.pickle')
        coastal_vulnerability.interpolate_sealevelrise_points(
            shore_point_path, slr_point_vector_path,
            slr_fieldname, target_pickle_path)

        with open(target_pickle_path, 'rb') as file:
            actual_values = pickle.load(file)
        expected_values = numpy.array([numpy.nan])

        numpy.testing.assert_almost_equal(
            actual_values.values(), expected_values, decimal=4)

    def test_long_aggregate_radius(self):
        """CV: handle an unreasonably long search radius in raster aggregation."""

        workspace_dir = self.workspace_dir
        raster_path = os.path.join(workspace_dir, 'simple_raster.tif')
        target_pickle_path = os.path.join(workspace_dir, 'target.pickle')
        sample_distance = 1.5

        # Make a simple raster
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N
        projection_wkt = srs.ExportToWkt()
        geotransform = [0, 0.5, 0.0, 0, 0.0, -0.5]
        n = 5
        nodata_val = -1
        gtiff_driver = gdal.GetDriverByName('GTiff')
        new_raster = gtiff_driver.Create(
            raster_path, n, n, 1, gdal.GDT_Int32, options=[
                'TILED=YES', 'BIGTIFF=YES', 'COMPRESS=LZW',
                'BLOCKXSIZE=16', 'BLOCKYSIZE=16'])
        new_raster.SetProjection(projection_wkt)
        new_raster.SetGeoTransform(geotransform)
        new_band = new_raster.GetRasterBand(1)
        array = numpy.array(range(n**2)).reshape((n, n))
        new_band.WriteArray(array)
        if nodata_val is not None:
            new_band.SetNoDataValue(nodata_val)
        new_raster.FlushCache()
        new_band = None
        new_raster = None

        # Make a vector proximate to the raster
        simple_points_path = os.path.join(workspace_dir, 'simple_points.shp')
        sampledata.create_vector_on_disk(
            [Point(0.1, -0.1),  # pixel (0,0): kernel origin out of bounds
             Point(1.25, -1.25),  # pixel (2,2): kernel origin & extent out of bounds
             Point(2.1, -2.1),  # pixel (4,4): kernel extent out of bounds
             ],
            srs.ExportToWkt(), filename=simple_points_path)

        coastal_vulnerability._aggregate_raster_values_in_radius(
            simple_points_path, raster_path, sample_distance,
            target_pickle_path, 'mean')

        with open(target_pickle_path, 'rb') as file:
            actual_values = pickle.load(file)

        expected_values = numpy.array([6.5454, 12.0, 17.4545])

        numpy.testing.assert_almost_equal(
            actual_values.values(), expected_values, decimal=4)

    def test_complete_run(self):
        """CV: regression test for a complete run with all optional arguments."""

        args = CoastalVulnerabilityTests.generate_base_args(self.workspace_dir)
        # these optional args aren't included in base_args:
        args['geomorphology_vector_path'] = os.path.join(
            INPUT_DATA, 'geomorphology_few_ranks.shp')
        args['geomorphology_fill_value'] = 3
        args['population_raster_path'] = os.path.join(
            INPUT_DATA, 'population.tif')
        args['population_radius'] = 16000
        args['slr_vector_path'] = os.path.join(
            INPUT_DATA, 'sea_level_rise.gpkg')
        args['slr_field'] = 'Trend'
        coastal_vulnerability.execute(args)

        actual_values_df = pandas.read_csv(
            os.path.join(args['workspace_dir'], 'coastal_exposure.csv'))
        expected_values_df = pandas.read_csv(
            os.path.join(REGRESSION_DATA, 'expected_coastal_exposure.csv'))
        pandas.testing.assert_frame_equal(actual_values_df, expected_values_df)

    def test_final_risk_calc(self):
        """CV: regression test for the final risk score calculation."""

        workspace_dir = self.workspace_dir
        base_point_vector_path = os.path.join(
            INPUT_DATA, "wwiii_shore_points_5000m.gpkg")

        risk_id_path_list = [
            (os.path.join(REGRESSION_DATA, 'expected_relief.pickle'), True, 'R_relief'),
            (os.path.join(REGRESSION_DATA, 'expected_wind.pickle'), True, 'R_wind'),
            (os.path.join(REGRESSION_DATA, 'expected_wave.pickle'), True, 'R_wave'),
            (os.path.join(REGRESSION_DATA, 'expected_surge.pickle'), True, None)
            ]
        habitat_protection_path = os.path.join(
            REGRESSION_DATA, 'expected_habitat_protection.csv')

        target_intermediate_vector_path = os.path.join(
            workspace_dir, 'intermediate_exposure.gpkg')
        target_point_vector_path = os.path.join(
            workspace_dir, 'coastal_exposure.gpkg')
        target_point_csv_path = os.path.join(
            workspace_dir, 'coastal_exposure.csv')

        coastal_vulnerability.assemble_results_and_calculate_exposure(
            risk_id_path_list, habitat_protection_path, base_point_vector_path,
            target_intermediate_vector_path, target_point_vector_path,
            target_point_csv_path)

        actual_values_df = pandas.read_csv(target_point_csv_path)
        expected_values_df = pandas.read_csv(
            os.path.join(REGRESSION_DATA, 'expected_final_risk.csv'))
        pandas.testing.assert_frame_equal(actual_values_df, expected_values_df)

    def test_geometric_mean_with_nan(self):
        """CV: test geometric mean function retuns `nan` with missing data."""
        array = numpy.array([1, 1, 1, 1, None], dtype=numpy.float)
        result = coastal_vulnerability._geometric_mean(array)
        assert numpy.isnan(result)

    def test_final_risk_calc_with_missing_data(self):
        """CV: test missing data at feature propogates to empty field in output."""
        target_vector_path = os.path.join(self.workspace_dir, 'target.gpkg')
        target_csv_path = os.path.join(self.workspace_dir, 'target.csv')
        # This gpkg has a feature with an empty field value for 'R_slr'
        shutil.copyfile(
            os.path.join(REGRESSION_DATA, 'test_missing_values.gpkg'),
            target_vector_path)

        coastal_vulnerability.calculate_final_risk(
            target_vector_path, target_csv_path)
        actual_values_df = pandas.read_csv(target_csv_path)

        # These fields should have missing values after the final calculations
        na_cols = ['exposure', 'habitat_role', 'exposure_no_habitats']
        na_data = [numpy.nan] * 3
        expected_df = pandas.DataFrame([na_data], columns=na_cols)

        pandas.testing.assert_frame_equal(
            actual_values_df[na_cols], expected_df)

    def test_binning_with_missing_data(self):
        """CV: test binning continuous values to ranks, w/ missing values."""
        n = 50
        mask = [10, 20, 30, 40]
        keys = range(n)
        numpy.random.seed(0)
        values = numpy.random.uniform(0, 1, n)
        values[mask] = numpy.nan
        missing_values_dict = dict(zip(keys, values))

        ranks_dict = coastal_vulnerability._bin_values_to_percentiles(
            missing_values_dict)

        # with random uniform values, all 5 ranks should be present
        expected_ranks = set([1, 2, 3, 4, 5])
        assert expected_ranks.issubset(ranks_dict.values())
        # and the masked indices should be nans
        assert all(numpy.isnan(numpy.array(ranks_dict.values())[mask]))

    def test_nodata_raster_aggregation(self):
        """CV: test raster aggregation over entirely nodata returns nan."""
        workspace_dir = self.workspace_dir
        raster_path = os.path.join(workspace_dir, 'nodata_raster.tif')
        target_pickle_path = os.path.join(workspace_dir, 'target.pickle')
        sample_distance = 1.5

        # Make a simple raster filled with all nodata
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N
        projection_wkt = srs.ExportToWkt()
        geotransform = [0, 0.5, 0.0, 0, 0.0, -0.5]
        n = 5
        nodata_val = -1
        gtiff_driver = gdal.GetDriverByName('GTiff')
        new_raster = gtiff_driver.Create(
            raster_path, n, n, 1, gdal.GDT_Int32, options=[
                'TILED=YES', 'BIGTIFF=YES', 'COMPRESS=LZW',
                'BLOCKXSIZE=16', 'BLOCKYSIZE=16'])
        new_raster.SetProjection(projection_wkt)
        new_raster.SetGeoTransform(geotransform)
        new_band = new_raster.GetRasterBand(1)
        array = numpy.array([nodata_val] * n**2).reshape((n, n))
        new_band.WriteArray(array)
        if nodata_val is not None:
            new_band.SetNoDataValue(nodata_val)
        new_raster.FlushCache()
        new_band = None
        new_raster = None

        # Make a vector proximate to the raster
        simple_points_path = os.path.join(workspace_dir, 'simple_points.shp')
        sampledata.create_vector_on_disk(
            [Point(0., 0.)],
            srs.ExportToWkt(), filename=simple_points_path)

        coastal_vulnerability._aggregate_raster_values_in_radius(
            simple_points_path, raster_path, sample_distance,
            target_pickle_path, 'density')

        with open(target_pickle_path, 'rb') as file:
            actual_values = pickle.load(file)

        expected_values = numpy.array([numpy.nan])
        numpy.testing.assert_almost_equal(
            actual_values.values(), expected_values, decimal=4)

    def test_exception_from_validate_polyline(self):
        """CV: raise ValueError on incorrect geometry type during validation.

        shelf_contour_vector_path must be a line geometry, here it's a polygon.
        """
        gpkg_driver = ogr.GetDriverByName("GPKG")
        shelf_poly_path = os.path.join(self.workspace_dir, 'shelf_poly.gpkg')
        vector = gpkg_driver.CreateDataSource(shelf_poly_path)
        vector.CreateLayer(
            'layer', osr.SpatialReference(), ogr.wkbPolygon)
        vector = None
        args = CoastalVulnerabilityTests.generate_base_args(self.workspace_dir)
        args['shelf_contour_vector_path'] = shelf_poly_path
        with self.assertRaises(ValueError):
            err_list = coastal_vulnerability.validate(args)
            for keys, err_strings in err_list:
                if 'Must be a polyline vector' in err_strings:
                    raise ValueError(err_list)

    def test_shore_points_on_single_polygon(self):
        """CV: test shore point creation with single polygon landmass."""
        args = CoastalVulnerabilityTests.generate_base_args(self.workspace_dir)
        aoi_path = os.path.join(self.workspace_dir, 'aoi.geojson')
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N
        wkt = srs.ExportToWkt()
        sampledata.create_vector_on_disk(
            [Polygon([(-200, -200), (200, -200), (200, 200), (-200, 200), (-200, -200)])],
            wkt, filename=aoi_path)

        landmass_path = os.path.join(
            self.workspace_dir, 'landmass.geojson')
        sampledata.create_vector_on_disk(
            [Polygon([(-100, -100), (100, -100), (100, 100), (-100, 100), (-100, -100)])],
            wkt, filename=landmass_path)

        args['aoi_vector_path'] = aoi_path
        args['landmass_vector_path'] = landmass_path
        args['model_resolution'] = 100

        polygon_pickle = os.path.join(
            self.workspace_dir, 'polygon.pickle')
        lines_pickle = os.path.join(
            self.workspace_dir, 'lines.pickle')
        lines_rtree = os.path.join(
            self.workspace_dir, 'rtree.dat')
        target_vector_path = os.path.join(
            self.workspace_dir, 'shore_points.gpkg')

        coastal_vulnerability.prepare_landmass_line_index(
            args['landmass_vector_path'], polygon_pickle,
            lines_pickle, lines_rtree)
        coastal_vulnerability.interpolate_shore_points(
            args['aoi_vector_path'], lines_pickle,
            args['model_resolution'], target_vector_path)

        vector = gdal.OpenEx(
            target_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
        layer = vector.GetLayer()
        n_points = layer.GetFeatureCount()
        assert n_points == 8

    def test_shore_points_on_multi_polygon(self):
        """CV: test shore point creation with multipolygon landmass."""

        args = CoastalVulnerabilityTests.generate_base_args(self.workspace_dir)
        aoi_path = os.path.join(self.workspace_dir, 'aoi.geojson')
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N
        wkt = srs.ExportToWkt()
        sampledata.create_vector_on_disk(
            [Polygon([(-200, -200), (200, -200), (200, 200), (-200, 200), (-200, -200)])],
            wkt, filename=aoi_path)

        landmass_path = os.path.join(
            self.workspace_dir, 'landmass.geojson')
        # two polygons disjoint from each other:
        sampledata.create_vector_on_disk(
            [Polygon([(-200, -200), (-100, -200), (-100, -100), (-200, -100), (-200, -200)]),
             Polygon([(100, 100), (200, 100), (200, 200), (100, 200), (100, 100)])],
            wkt, filename=landmass_path)

        args['aoi_vector_path'] = aoi_path
        args['landmass_vector_path'] = landmass_path
        args['model_resolution'] = 100

        polygon_pickle = os.path.join(
            self.workspace_dir, 'polygon.pickle')
        lines_pickle = os.path.join(
            self.workspace_dir, 'lines.pickle')
        lines_rtree = os.path.join(
            self.workspace_dir, 'rtree.dat')
        target_vector_path = os.path.join(
            self.workspace_dir, 'shore_points.gpkg')

        coastal_vulnerability.prepare_landmass_line_index(
            args['landmass_vector_path'], polygon_pickle,
            lines_pickle, lines_rtree)
        coastal_vulnerability.interpolate_shore_points(
            args['aoi_vector_path'], lines_pickle,
            args['model_resolution'], target_vector_path)

        vector = gdal.OpenEx(
            target_vector_path, gdal.OF_VECTOR | gdal.GA_ReadOnly)
        layer = vector.GetLayer()
        n_points = layer.GetFeatureCount()
        assert n_points == 8

    def test_no_wwiii_coverage(self):
        """CV: test exception when shore points are outside max wwiii distance."""
        args = CoastalVulnerabilityTests.generate_base_args(self.workspace_dir)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(26910)  # UTM Zone 10N

        simple_points_path = os.path.join(self.workspace_dir, 'simple_points.shp')
        sampledata.create_vector_on_disk(
            [Point(0.0, 0.0)],
            srs.ExportToWkt(), filename=simple_points_path)

        with self.assertRaises(ValueError):
            coastal_vulnerability.interpolate_wwiii_to_shore(
                simple_points_path, args['wwiii_vector_path'],
                'target.gpkg')


def assert_pickled_arrays_almost_equal(
        actual_values_pickle_path, expected_values_pickle_path):
    """Open two pickled dicts and assert keys and values match."""
    with open(actual_values_pickle_path, 'rb') as pickle_file:
        actual_values_dict = pickle.load(pickle_file)
    actual_values = actual_values_dict.values()
    actual_fids = actual_values_dict.keys()

    with open(expected_values_pickle_path, 'rb') as pickle_file:
        expected_values_dict = pickle.load(pickle_file)
    expected_values = expected_values_dict.values()
    expected_fids = expected_values_dict.keys()

    numpy.testing.assert_array_almost_equal(
        actual_values, expected_values, decimal=2)
    numpy.testing.assert_array_equal(actual_fids, expected_fids)
