def test_rectangular(self):
        lons = numpy.array(range(100)).reshape((10, 10))
        lats = numpy.negative(lons)

        mesh = RectangularMesh(lons, lats, depths=None)
        bounding_mesh = mesh._get_bounding_mesh()
        expected_lons = numpy.array([
            0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
            19, 29, 39, 49, 59, 69, 79, 89,
            99, 98, 97, 96, 95, 94, 93, 92, 91,
            90, 80, 70, 60, 50, 40, 30, 20, 10
        ])
        expected_lats = numpy.negative(expected_lons)
        self.assertTrue((bounding_mesh.lons == expected_lons).all())
        self.assertTrue((bounding_mesh.lats == expected_lats).all())
        self.assertIsNone(bounding_mesh.depths)

        depths = lons + 10
        mesh = RectangularMesh(lons, lats, depths)
        expected_depths = expected_lons + 10
        bounding_mesh = mesh._get_bounding_mesh()
        self.assertIsNotNone(bounding_mesh.depths)
        self.assertTrue((bounding_mesh.depths
                         == expected_depths.flatten()).all())

        bounding_mesh = mesh._get_bounding_mesh(with_depths=False)
        self.assertIsNone(bounding_mesh.depths)

    def test_simple(self):
        lons = numpy.array([numpy.arange(-1, 1.2, 0.2)] * 11)
        lats = lons.transpose() + 1
        depths = lats + 10
        mesh = RectangularMesh(lons, lats, depths)

        check = lambda lon, lat, depth, expected_distance, **kwargs: \
            self.assertAlmostEqual(
                mesh.get_joyner_boore_distance(
                    Mesh.from_points_list([Point(lon, lat, depth)])
                )[0],
                expected_distance, **kwargs
            )

        check(lon=0, lat=0.5, depth=0, expected_distance=0)
        check(lon=1, lat=1, depth=0, expected_distance=0)
        check(lon=0.6, lat=-1, depth=0,
              expected_distance=Point(0.6, -1).distance(Point(0.6, 0)),
              delta=0.1)
        check(lon=-0.8, lat=2.1, depth=10,
              expected_distance=Point(-0.8, 2.1).distance(Point(-0.8, 2)),
              delta=0.02)
        check(lon=0.75, lat=2.3, depth=3,
              expected_distance=Point(0.75, 2.3).distance(Point(0.75, 2)),
              delta=0.04)

 def test_even_rows_even_columns_with_depths(self):
     lons = numpy.array([[10, 20], [12, 22]])
     lats = numpy.array([[10, -10], [8, -9]])
     depths = numpy.array([[2, 3], [4, 5]])
     mesh = RectangularMesh(lons, lats, depths=depths)
     self.assertEqual(mesh.get_middle_point(),
                      Point(15.996712, -0.250993, 3.5))

 def test_mesh_of_one_point(self):
     lons = numpy.array([[1]])
     lats = numpy.array([[0]])
     depths = numpy.array([[1]])
     mesh = RectangularMesh(lons, lats, depths)
     target_mesh = Mesh.from_points_list([Point(1, 0), Point(0.5, 0)])
     dists = mesh.get_joyner_boore_distance(target_mesh)
     expected_dists = [0, Point(0.5, 0).distance(Point(1, 0))]
     self.assertTrue(numpy.allclose(dists, expected_dists, atol=0.2))

 def _test(self, lons, lats, depths, expected_coords):
     mesh = RectangularMesh(lons, lats, depths)
     proj, polygon = mesh._get_proj_enclosing_polygon()
     self.assertTrue(polygon.is_valid)
     self.assertEqual(list(polygon.interiors), [])
     coords = numpy.array(proj(*numpy.array(polygon.exterior).transpose(),
                               reverse=True)).transpose()
     numpy.testing.assert_almost_equal(coords, expected_coords, decimal=4)
     return polygon

 def test_mesh_of_two_points(self):
     lons = numpy.array([[0, 0.5, 1]], float)
     lats = numpy.array([[0, 0, 0]], float)
     depths = numpy.array([[1, 0, 1]], float)
     mesh = RectangularMesh(lons, lats, depths)
     target_mesh = Mesh.from_points_list([Point(0.5, 1), Point(0.5, 0)])
     dists = mesh.get_joyner_boore_distance(target_mesh)
     expected_dists = [Point(0.5, 1).distance(Point(0.5, 0)), 0]
     numpy.testing.assert_almost_equal(dists, expected_dists)

 def _test(self, points, site, expected_distance):
     lons, lats, depths = numpy.array(points).transpose()
     lons = lons.transpose()
     lats = lats.transpose()
     depths = depths.transpose()
     mesh = RectangularMesh(lons, lats, depths)
     distance = mesh.get_joyner_boore_distance(
         Mesh.from_points_list([Point(*site)])
     )[0]
     self.assertAlmostEqual(distance, expected_distance, delta=0.02)

 def test_mesh_width(self):
     lons = numpy.array([[0.1, 0.1, 0.1, 0.1],
                         [0.1, 0.1, 0.1, 0.1],
                         [0.1, 0.1, 0.1, 0.1]])
     lats = numpy.array([[0.1, 0.10899322, 0.11798643, 0.12697965],
                         [0.1, 0.10899322, 0.11798643, 0.12697965],
                         [0.1, 0.10899322, 0.11798643, 0.12697965]])
     depths = numpy.array([[2.0, 2.0, 2.0, 2.0],
                           [3.0, 3.0, 3.0, 3.0],
                           [4.0, 4.0, 4.0, 4.0]])
     mesh = RectangularMesh(lons, lats, depths)
     self.assertAlmostEqual(mesh.get_mean_width(), 2.0)

 def test_vertical_mesh(self):
     lons = numpy.array([[0, 1, 2], [0, 1, 2]])
     lats = numpy.array([[0, 0, 0], [0, 0, 0]])
     depths = numpy.array([[1, 1, 1], [2, 2, 2]])
     mesh = RectangularMesh(lons, lats, depths)
     target_mesh = Mesh.from_points_list([Point(0.5, 0), Point(0.5, 1),
                                          Point(0.5, 5)])
     dists = mesh.get_joyner_boore_distance(target_mesh)
     expected_dists = [
         0, Point(0.5, 1).distance(Point(0.5, 0)),
         Point(0.5, 5).distance(Point(0.5, 0))
     ]
     self.assertTrue(numpy.allclose(dists, expected_dists, atol=3))

    def test_single_column(self):
        lons = numpy.array([[0], [1], [2], [3], [4], [5]])
        lats = numpy.array([[-1], [-2], [-3], [-4], [-5], [-6]])
        mesh = RectangularMesh(lons, lats, depths=None)
        bounding_mesh = mesh._get_bounding_mesh()
        self.assertTrue((bounding_mesh.lons == lons.flatten()).all())
        self.assertTrue((bounding_mesh.lats == lats.flatten()).all())
        self.assertIsNone(bounding_mesh.depths)

        depths = numpy.array([[10], [11], [12], [13], [14], [15]])
        mesh = RectangularMesh(lons, lats, depths)
        bounding_mesh = mesh._get_bounding_mesh()
        self.assertIsNotNone(bounding_mesh.depths)
        self.assertTrue((bounding_mesh.depths == depths.flatten()).all())

        bounding_mesh = mesh._get_bounding_mesh(with_depths=False)
        self.assertIsNone(bounding_mesh.depths)

    def test_single_row(self):
        lons = numpy.array([[0, 1, 2, 3, 4, 5]])
        lats = numpy.array([[-1, -2, -3, -4, -5, -6]])
        mesh = RectangularMesh(lons, lats, depths=None)
        bounding_mesh = mesh._get_bounding_mesh()
        self.assertIsInstance(bounding_mesh, Mesh)
        self.assertTrue((bounding_mesh.lons == lons[0]).all())
        self.assertTrue((bounding_mesh.lats == lats[0]).all())
        self.assertIsNone(bounding_mesh.depths)

        depths = numpy.array([[10, 11, 12, 13, 14, 15]])
        mesh = RectangularMesh(lons, lats, depths)
        bounding_mesh = mesh._get_bounding_mesh()
        self.assertIsNotNone(bounding_mesh.depths)
        self.assertTrue((bounding_mesh.depths == depths[0]).all())

        bounding_mesh = mesh._get_bounding_mesh(with_depths=False)
        self.assertIsNone(bounding_mesh.depths)

    def test_dip_over_90_degree(self):
        top = [Point(0, -0.01), Point(0, 0.01)]
        bottom = [Point(-0.01, -0.01, 1.11), Point(-0.01, 0.01, 1.11)]

        mesh = RectangularMesh.from_points_list([top, bottom])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        # dip must be still in a range 0..90
        self.assertAlmostEqual(dip, 45, delta=0.05)
        # strike must be reversed
        self.assertAlmostEqual(strike, 180, delta=0.05)

    def from_fault_data(cls, fault_trace, upper_seismogenic_depth,
                        lower_seismogenic_depth, dip, mesh_spacing):
        """
        Create and return a fault surface using fault source data.

        :param fault_trace:
            Geographical line representing the intersection between
            the fault surface and the earth surface, an instance
            of :class:`nhlib.Line`.
        :param upper_seismo_depth:
            Minimum depth ruptures can reach, in km (i.e. depth
            to fault's top edge).
        :param lower_seismo_depth:
            Maximum depth ruptures can reach, in km (i.e. depth
            to fault's bottom edge).
        :param dip:
            Dip angle (i.e. angle between fault surface
            and earth surface), in degrees.
        :param mesh_spacing:
            Distance between two subsequent points in a mesh, in km.
        :returns:
            An instance of :class:`SimpleFaultSurface` created using that data.

        Uses :meth:`check_fault_data` for checking parameters.
        """
        cls.check_fault_data(fault_trace, upper_seismogenic_depth,
                             lower_seismogenic_depth, dip, mesh_spacing)
        # Loops over points in the top edge, for each point
        # on the top edge compute corresponding point on the bottom edge, then
        # computes equally spaced points between top and bottom points.

        vdist_top = upper_seismogenic_depth
        vdist_bottom = lower_seismogenic_depth

        hdist_top = vdist_top / math.tan(math.radians(dip))
        hdist_bottom = vdist_bottom / math.tan(math.radians(dip))

        strike = fault_trace[0].azimuth(fault_trace[-1])
        azimuth = (strike + 90.0) % 360

        mesh = []
        for point in fault_trace.resample(mesh_spacing):
            top = point.point_at(hdist_top, vdist_top, azimuth)
            bottom = point.point_at(hdist_bottom, vdist_bottom, azimuth)
            mesh.append(top.equally_spaced_points(bottom, mesh_spacing))

        # number of rows corresponds to number of points along dip
        # number of columns corresponds to number of points along strike
        surface_points = numpy.array(mesh).transpose().tolist()
        mesh = RectangularMesh.from_points_list(surface_points)
        assert 1 not in mesh.shape
        return cls(mesh)

    def test_one_cell(self):
        top = [Point(0, -0.01), Point(0, 0.01)]
        bottom = [Point(0.01, -0.01, 1.11), Point(0.01, 0.01, 1.11)]

        mesh = RectangularMesh.from_points_list([top, bottom])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, 45, delta=0.05)
        self.assertAlmostEqual(strike, 0, delta=0.05)

        row1 = [Point(45, -0.1), Point(45.2, 0.1)]
        row2 = [Point(45, -0.1, 1), Point(45.2, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, 90)
        self.assertAlmostEqual(strike, 45, delta=0.1)

        row1 = [Point(90, -0.1), Point(90, 0.1)]
        row2 = [Point(90, -0.1, 1), Point(90, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, 90)
        self.assertAlmostEqual(strike, 0, delta=0.1)

    def test_one_cell_unequal_area(self):
        # top-left triangle is vertical, has dip of 90 degrees, zero
        # strike and area of 1 by 1 over 2. bottom-right one has dip
        # of atan2(1, sqrt(2) / 2.0) which is 54.73561 degrees, strike
        # of 45 degrees and area that is 1.73246136 times area of the
        # first one's. weighted mean dip is 67.5 degrees and weighted
        # mean strike is 28.84 degrees
        top = [Point(0, -0.01), Point(0, 0.01)]
        bottom = [Point(0, -0.01, 2.22), Point(0.02, 0.01, 2.22)]

        mesh = RectangularMesh.from_points_list([top, bottom])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, 67.5, delta=0.05)
        self.assertAlmostEqual(strike, 28.84, delta=0.05)

    def test_from_points_list(self):
        lons = [[0, 1], [2, 3], [4, 5]]
        lats = [[1, 2], [-1, -2], [10, 20]]
        depths = [[11.1, 11.2], [11.3, 11.4], [11.5, 11.6]]
        points = [
            [Point(lons[i][j], lats[i][j], depths[i][j])
             for j in xrange(len(lons[i]))]
            for i in xrange(len(lons))
        ]
        mesh = RectangularMesh.from_points_list(points)
        self.assertTrue((mesh.lons == lons).all())
        self.assertTrue((mesh.lats == lats).all())
        self.assertTrue((mesh.depths == depths).all())

        points = [
            [Point(lons[i][j], lats[i][j], depth=0)
             for j in xrange(len(lons[i]))]
            for i in xrange(len(lons))
        ]
        mesh = RectangularMesh.from_points_list(points)
        self.assertTrue((mesh.lons == lons).all())
        self.assertTrue((mesh.lats == lats).all())
        self.assertIsNone(mesh.depths)

    def from_fault_data(cls, edges, mesh_spacing):
        """
        Create and return a fault surface using fault source data.

        :param edges:
            A list of at least two horizontal edges of the surface
            as instances of :class:`nhlib.geo.line.Line`. The list
            should be in top-to-bottom order (the shallowest edge
            first).
        :param mesh_spacing:
            Distance between two subsequent points in a mesh, in km.
        :returns:
            An instance of :class:`ComplexFaultSurface` created using
            that data.
        :raises ValueError:
            If requested mesh spacing is too big for the surface geometry
            (doesn't allow to put a single mesh cell along length and/or
            width).

        Uses :meth:`check_fault_data` for checking parameters.
        """
        cls.check_fault_data(edges, mesh_spacing)

        mean_length = numpy.mean([edge.get_length() for edge in edges])
        num_hor_points = int(round(mean_length / mesh_spacing)) + 1
        if num_hor_points <= 1:
            raise ValueError(
                'mesh spacing %.1f km is too big for mean length %.1f km' %
                (mesh_spacing, mean_length)
            )
        edges = [edge.resample_to_num_points(num_hor_points).points
                 for i, edge in enumerate(edges)]

        vert_edges = [Line(v_edge) for v_edge in zip(*edges)]
        mean_width = numpy.mean([v_edge.get_length() for v_edge in vert_edges])
        num_vert_points = int(round(mean_width / mesh_spacing)) + 1
        if num_vert_points <= 1:
            raise ValueError(
                'mesh spacing %.1f km is too big for mean width %.1f km' %
                (mesh_spacing, mean_width)
            )

        points = zip(*[v_edge.resample_to_num_points(num_vert_points).points
                       for v_edge in vert_edges])
        mesh = RectangularMesh.from_points_list(points)
        assert 1 not in mesh.shape
        return cls(mesh)

    def test_on_surface(self):
        row1 = [Point(0, 0), Point(0, 1)]
        row2 = [Point(1, 0), Point(1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 0, delta=0.5)

        row1 = [Point(0, 0), Point(0, -1)]
        row2 = [Point(1, 0), Point(1, -1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 180, delta=0.5)

        row1 = [Point(0, 0), Point(1, 1)]
        row2 = [Point(1, 0), Point(2, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 45, delta=0.01)

        row1 = [Point(0, 0), Point(1, -1)]
        row2 = [Point(1, 0), Point(2, -1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 135, delta=0.01)

        row1 = [Point(0, 0), Point(-1, -1)]
        row2 = [Point(-1, 0), Point(-2, -1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 225, delta=0.01)

        row1 = [Point(0, 0), Point(-1, 1)]
        row2 = [Point(-1, 0), Point(-2, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertEqual(dip, 0)
        self.assertAlmostEqual(strike, 315, delta=0.01)

    def test_two_cells(self):
        top = [Point(0, -0.01), Point(0, 0.01)]
        middle = [Point(0.01, -0.01, 1.11), Point(0.01, 0.01, 1.11)]
        bottom = [Point(0.01, -0.01, 2.22), Point(0.01, 0.01, 2.22)]

        mesh = RectangularMesh.from_points_list([top, middle, bottom])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, math.degrees(math.atan2(2, 1)), delta=0.1)
        self.assertAlmostEqual(strike, 0, delta=0.02)

        bottom = [Point(0.01, -0.01, 3.33), Point(0.01, 0.01, 3.33)]
        mesh = RectangularMesh.from_points_list([top, middle, bottom])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, math.degrees(math.atan2(3, 1)), delta=0.1)
        self.assertAlmostEqual(strike, 0, delta=0.02)

        row1 = [Point(90, -0.1), Point(90, 0), Point(90, 0.1)]
        row2 = [Point(90, -0.1, 1), Point(90, 0, 1), Point(90, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(dip, 90)
        assert_angles_equal(self, strike, 360, delta=1e-7)

        row1 = [Point(-90.1, -0.1), Point(-90, 0), Point(-89.9, 0.1)]
        row2 = [Point(-90.0, -0.1, 1), Point(-89.9, 0, 1), Point(-89.8, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(strike, 45, delta=1e-4)

        row1 = [Point(-90.1, -0.1), Point(-90, 0), Point(-89.9, 0.1)]
        row2 = [Point(-90.0, -0.1, 1), Point(-89.9, 0, 1), Point(-89.8, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(strike, 45, delta=1e-3)

        row1 = [Point(-90.1, -0.1), Point(-90, 0), Point(-89.9, 0.1)]
        row2 = [Point(-90.2, -0.1, 1), Point(-90.1, 0, 1), Point(-90, 0.1, 1)]
        mesh = RectangularMesh.from_points_list([row1, row2])
        dip, strike = mesh.get_mean_inclination_and_azimuth()
        self.assertAlmostEqual(strike, 225, delta=1e-3)

 def test_odd_rows_odd_columns_with_depths(self):
     lons = numpy.array([numpy.arange(-1, 1.2, 0.2)] * 11)
     lats = lons.transpose() * 10
     depths = lats + 10
     mesh = RectangularMesh(lons, lats, depths)
     self.assertEqual(mesh.get_middle_point(), Point(0, 0, 10))