Python point.Point类代码示例

    def circular_distance_from_point(self, point, distance, **kwargs):
        '''
        Select earthquakes within a distance from a Point

        :param point:
            Centre point as instance of nhlib.geo.point.Point class

        :param float distance:
            Distance (km)

        :returns:
            Instance of :class:`openquake.hmtk.seismicity.catalogue.Catalogue`
            containing only selected events
        '''

        if kwargs['distance_type'] is 'epicentral':
            locations = Mesh(
                self.catalogue.data['longitude'],
                self.catalogue.data['latitude'],
                np.zeros(len(self.catalogue.data['longitude']), dtype=float))
            point = Point(point.longitude, point.latitude, 0.0)
        else:
            locations = self.catalogue.hypocentres_as_mesh()

        is_close = point.closer_than(locations, distance)

        return self.select_catalogue(is_close)

    def cartesian_square_centred_on_point(self, point, distance, **kwargs):
        '''
        Select earthquakes from within a square centered on a point

        :param point:
            Centre point as instance of nhlib.geo.point.Point class

        :param distance:
            Distance (km)

        :returns:
            Instance of :class:`openquake.hmtk.seismicity.catalogue.Catalogue`
            class containing only selected events
        '''
        point_surface = Point(point.longitude, point.latitude, 0.)
        # As distance is
        north_point = point_surface.point_at(distance, 0., 0.)
        east_point = point_surface.point_at(distance, 0., 90.)
        south_point = point_surface.point_at(distance, 0., 180.)
        west_point = point_surface.point_at(distance, 0., 270.)
        is_long = np.logical_and(
            self.catalogue.data['longitude'] >= west_point.longitude,
            self.catalogue.data['longitude'] < east_point.longitude)
        is_surface = np.logical_and(
            is_long,
            self.catalogue.data['latitude'] >= south_point.latitude,
            self.catalogue.data['latitude'] < north_point.latitude)

        upper_depth, lower_depth = _check_depth_limits(kwargs)
        is_valid = np.logical_and(
            is_surface,
            self.catalogue.data['depth'] >= upper_depth,
            self.catalogue.data['depth'] < lower_depth)

        return self.select_catalogue(is_valid)

 def setUp(self):
     '''
     '''
     self.fault = None
     self.regionalisation = None
     self.msr = [(WC1994(), 1.0)]
     self.msr_sigma = [(-1.5, 0.15), (0.0, 0.7), (1.5, 0.15)]
     self.shear_mod = [(30.0, 0.8), (35.0, 0.2)]
     self.dlr = [(1.25E-5, 1.0)]
     self.config = [{}]
     self.slip = [(10.0, 1.0)]
     x0 = Point(30., 30., 0.)
     x1 = x0.point_at(30., 0., 30.)
     x2 = x1.point_at(30., 0., 60.)
     # Total length is 60 km
     self.trace = Line([x0, x1, x2])
     self.dip = 90.
     self.upper_depth = 0.
     self.lower_depth = 20.
     self.simple_fault = SimpleFaultGeometry(self.trace,
                                             self.dip,
                                             self.upper_depth,
                                             self.lower_depth)
             # Creates a trace ~60 km long made of 3 points
     upper_edge = Line([x0, x1, x2])
     lower_edge = Line([x0.point_at(40., 20., 130.),
                        x1.point_at(42., 25., 130.),
                        x2.point_at(41., 22., 130.)])
     self.complex_fault = ComplexFaultGeometry([upper_edge, lower_edge],
                                               2.0)

def build_planar_surface(geometry):
    """
    Builds the planar rupture surface from the openquake.nrmllib.models
    instance
    """
    # Read geometry from wkt
    geom = wkt.loads(geometry.wkt)
    top_left = Point(geom.xy[0][0],
                     geom.xy[1][0],
                     geometry.upper_seismo_depth)
    top_right = Point(geom.xy[0][1],
                      geom.xy[1][1],
                      geometry.upper_seismo_depth)
    strike = top_left.azimuth(top_right)
    dip_dir = (strike + 90.) % 360.
    depth_diff = geometry.lower_seismo_depth - geometry.upper_seismo_depth
    bottom_right = top_right.point_at(
        depth_diff / np.tan(geometry.dip * (np.pi / 180.)),
        depth_diff,
        dip_dir)
    bottom_left = top_left.point_at(
        depth_diff / np.tan(geometry.dip * (np.pi / 180.)),
        depth_diff,
        dip_dir)
    return PlanarSurface(1.0,
                         strike,
                         geometry.dip,
                         top_left,
                         top_right,
                         bottom_right,
                         bottom_left)

def _rup_to_point(distance, surface, origin, azimuth, distance_type='rjb',
        iter_stop=1E-5, maxiter=1000):
    """

    """
    pt0 = origin
    pt1 = origin.point_at(distance, 0., azimuth)
    r_diff = np.inf
    iterval = 0
    while (np.fabs(r_diff) >= iter_stop) and (iterval <= maxiter):
        pt1mesh = Mesh(np.array([pt1.longitude]),
                       np.array([pt1.latitude]),
                       None)
        if distance_type == 'rjb':
            r_diff =  distance - surface.get_joyner_boore_distance(pt1mesh)
        elif distance_type == 'rrup':
            r_diff =  distance - surface.get_min_distance(pt1mesh)
        else:
            raise ValueError('Distance type must be rrup or rjb!')
        pt0 = Point(pt1.longitude, pt1.latitude)
        if r_diff > 0.:
            pt1 = pt0.point_at(r_diff, 0., azimuth)
        else:
            pt1 = pt0.point_at(r_diff, 0., (azimuth + 180.) % 360.)
    return pt1

 def setUp(self):
     '''
     Creates a complex fault typology
     '''
     x0 = Point(30., 30., 0.)
     x1 = x0.point_at(30., 0., 30.)
     x2 = x1.point_at(30., 0., 60.)
     upper_edge = Line([x0, x1, x2])
     lower_edge = Line([x0.point_at(40., 20., 130.),
                        x1.point_at(42., 25., 130.),
                        x2.point_at(41., 22., 130.)])
     self.edges = [upper_edge, lower_edge]
     self.fault = None

    def check_surface_validity(cls, edges):
        """
        Check validity of the surface.

        Project edge points to vertical plane anchored to surface upper left
        edge and with strike equal to top edge strike. Check that resulting
        polygon is valid.

        This method doesn't have to be called by hands before creating the
        surface object, because it is called from :meth:`from_fault_data`.
        """
        # extract coordinates of surface boundary (as defined from edges)
        full_boundary = []
        left_boundary = []
        right_boundary = []

        for i in range(1, len(edges) - 1):
            left_boundary.append(edges[i].points[0])
            right_boundary.append(edges[i].points[-1])

        full_boundary.extend(edges[0].points)
        full_boundary.extend(right_boundary)
        full_boundary.extend(edges[-1].points[::-1])
        full_boundary.extend(left_boundary[::-1])

        lons = [p.longitude for p in full_boundary]
        lats = [p.latitude for p in full_boundary]
        depths = [p.depth for p in full_boundary]

        # define reference plane. Corner points are separated by an arbitrary
        # distance of 10 km. The mesh spacing is set to 2 km. Both corner
        # distance and mesh spacing values do not affect the algorithm results.
        ul = edges[0].points[0]
        strike = ul.azimuth(edges[0].points[-1])
        dist = 10.
        mesh_spacing = 2.

        ur = ul.point_at(dist, 0, strike)
        bl = Point(ul.longitude, ul.latitude, ul.depth + dist)
        br = bl.point_at(dist, 0, strike)

        # project surface boundary to reference plane and check for
        # validity.
        ref_plane = PlanarSurface.from_corner_points(
            mesh_spacing, ul, ur, br, bl
        )
        _, xx, yy = ref_plane._project(lons, lats, depths)
        coords = [(x, y) for x, y in zip(xx, yy)]
        p = shapely.geometry.Polygon(coords)
        if not p.is_valid:
            raise ValueError('Edges points are not in the right order')

def assert_mesh_is(testcase, surface, expected_mesh):
    mesh = surface.get_mesh()
    testcase.assertIs(mesh, surface.get_mesh())

    expected_mesh = list(itertools.chain(*expected_mesh))
    testcase.assertEqual(len(mesh), len(expected_mesh))
    testcase.assertIsInstance(mesh, Mesh)

    for i, point in enumerate(mesh):
        expected_point = Point(*expected_mesh[i])
        distance = expected_point.distance(point) * 1e3

        testcase.assertAlmostEqual(
            0, distance, delta=2, # allow discrepancy of 2 meters
            msg="point %d is off: %s != %s (distance is %.3fm)"
                % (i, point, expected_point, distance)
        )

    def getLength(self):
        """
        Compute length of rupture (km). For EdgeRupture, we compute the length
        as the length of the top edge projected to the surface.

        Returns:
            float: Rupture length in km.
        """
        lons = self._toplons
        lats = self._toplats
        seg = self._group_index
        groups = np.unique(seg)
        ng = len(groups)
        rlength = 0
        for i in range(ng):
            group_segments = np.where(groups[i] == seg)[0]
            nseg = len(group_segments) - 1
            for j in range(nseg):
                ind = group_segments[j]
                P0 = Point(lons[ind], lats[ind])
                P1 = Point(lons[ind + 1], lats[ind + 1])
                dist = P0.distance(P1)
                rlength = rlength + dist
        return rlength

    def test_build_fault_model(self):
        # Tests the constuction of a fault model with two faults (1 simple,
        # 1 complex) each with two mfd rates - should produce four sources
        self.model = mtkActiveFaultModel('001', 'A Fault Model', faults=[])
        x0 = Point(30., 30., 0.)
        x1 = x0.point_at(30., 0., 30.)
        x2 = x1.point_at(30., 0., 60.)
        # Total length is 60 km
        trace = Line([x0, x1, x2])
        simple_fault = SimpleFaultGeometry(trace, 90., 0., 20.)
        # Creates a trace ~60 km long made of 3 points
        upper_edge = Line([x0, x1, x2])
        lower_edge = Line(
            [x0.point_at(40., 20., 130.),
             x1.point_at(42., 25., 130.),
             x2.point_at(41., 22., 130.)])
        complex_fault = ComplexFaultGeometry([upper_edge, lower_edge], 2.0)
        config = [{'MFD_spacing': 0.1,
                   'Maximum_Magnitude': 7.0,
                   'Maximum_Uncertainty': None,
                   'Model_Name': 'Characteristic',
                   'Model_Weight': 0.5,
                   'Sigma': 0.1,
                   'Lower_Bound': -1.,
                   'Upper_Bound': 1.},
                  {'MFD_spacing': 0.1,
                   'Maximum_Magnitude': 7.5,
                   'Maximum_Uncertainty': None,
                   'Model_Name': 'Characteristic',
                   'Model_Weight': 0.5,
                   'Sigma': 0.1,
                   'Lower_Bound': -1.,
                   'Upper_Bound': 1.}]
        fault1 = mtkActiveFault('001', 'Simple Fault 1', simple_fault,
                                [(10.0, 1.0)], -90., None,
                                aspect_ratio=1.0,
                                scale_rel=[(WC1994(), 1.0)],
                                shear_modulus=[(30.0, 1.0)],
                                disp_length_ratio=[(1E-5, 1.0)])
        fault1.generate_config_set(config)
        fault2 = mtkActiveFault('002', 'Complex Fault 1', complex_fault,
                                [(10.0, 1.0)], -90., None,
                                aspect_ratio=1.0,
                                scale_rel=[(WC1994(), 1.0)],
                                shear_modulus=[(30.0, 1.0)],
                                disp_length_ratio=[(1E-5, 1.0)])
        fault2.generate_config_set(config)
        self.model.faults = [fault1, fault2]

        # Generate source model
        self.model.build_fault_model()
        self.assertEqual(len(self.model.source_model.sources), 4)
        # First source should be an instance of a mtkSimpleFaultSource
        model1 = self.model.source_model.sources[0]
        self.assertTrue(isinstance(model1, mtkSimpleFaultSource))
        self.assertEqual(model1.id, '001_1')
        self.assertAlmostEqual(model1.mfd.min_mag, 6.9)
        np.testing.assert_array_almost_equal(
            np.log10(np.array(model1.mfd.occurrence_rates)),
            np.array([-2.95320041, -2.54583708, -2.953200413]))

        # Second source should be an instance of a mtkSimpleFaultSource
        model2 = self.model.source_model.sources[1]
        self.assertTrue(isinstance(model2, mtkSimpleFaultSource))
        self.assertEqual(model2.id, '001_2')
        self.assertAlmostEqual(model2.mfd.min_mag, 7.4)
        np.testing.assert_array_almost_equal(
            np.log10(np.array(model2.mfd.occurrence_rates)),
            np.array([-3.70320041, -3.29583708, -3.70320041]))

        # Third source should be an instance of a mtkComplexFaultSource
        model3 = self.model.source_model.sources[2]
        self.assertTrue(isinstance(model3, mtkComplexFaultSource))
        self.assertEqual(model3.id, '002_1')
        self.assertAlmostEqual(model3.mfd.min_mag, 6.9)
        np.testing.assert_array_almost_equal(
            np.log10(np.array(model3.mfd.occurrence_rates)),
            np.array([-2.59033387, -2.18297054, -2.59033387]))

        # Fourth source should be an instance of a mtkComplexFaultSource
        model4 = self.model.source_model.sources[3]
        self.assertTrue(isinstance(model4, mtkComplexFaultSource))
        self.assertEqual(model4.id, '002_2')
        self.assertAlmostEqual(model4.mfd.min_mag, 7.4)
        np.testing.assert_array_almost_equal(
            np.log10(np.array(model4.mfd.occurrence_rates)),
            np.array([-3.34033387, -2.93297054, -3.34033387]))