Source code for openquake.hazardlib.geo.point

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2023 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.

"""
Module :mod:`openquake.hazardlib.geo.point` defines :class:`Point`.
"""
import numpy
import shapely.geometry

from openquake.hazardlib.geo import geodetic
from openquake.hazardlib.geo import utils as geo_utils


[docs]class Point(object): """ This class represents a geographical point in terms of longitude, latitude, and depth (with respect to the Earth surface). :param longitude: Point longitude, in decimal degrees. :type longitude: float :param latitude: Point latitude, in decimal degrees. :type latitude: float :param depth: Point depth (default to 0.0), in km. Depth > 0 indicates a point below the earth surface, and depth < 0 above the earth surface. :type depth: float """ #: The distance between two points for them to be considered equal, #: in km. EQUALITY_DISTANCE = 1e-3 def __init__(self, longitude, latitude, depth=0.0): if not depth < geodetic.EARTH_RADIUS: raise ValueError("The depth must be less than " "the Earth's radius (6371.0 km)") if not depth > geodetic.EARTH_ELEVATION: raise ValueError("The depth must be greater than the maximum " "elevation on Earth's surface (-8.848 km)") if not -180.0 <= longitude <= 180.0: raise ValueError("longitude %.6f outside range" % longitude) if not -90.0 <= latitude <= 90.0: raise ValueError("latitude %.6f outside range" % latitude) self.depth = float(depth) self.latitude = float(latitude) self.longitude = float(longitude) @property def x(self): """Alias for .longitude""" return self.longitude @property def y(self): """Alias for .latitude""" return self.latitude @property def z(self): """Alias for .depth""" return self.depth @property def wkt2d(self): """ Generate WKT (Well-Known Text) to represent this point in 2 dimensions (ignoring depth). """ return 'POINT(%s %s)' % (self.longitude, self.latitude)
[docs] def round(self, digits=5): """ :returns: a new Point with lon, lat, depth rounded to 5 digits """ lon = round(self.longitude, digits) lat = round(self.latitude, digits) dep = round(self.depth, digits) return Point(lon, lat, dep)
[docs] def point_at(self, horizontal_distance, vertical_increment, azimuth): """ Compute the point with given horizontal, vertical distances and azimuth from this point. :param horizontal_distance: Horizontal distance, in km. :type horizontal_distance: float :param vertical_increment: Vertical increment, in km. When positive, the new point has a greater depth. When negative, the new point has a smaller depth. :type vertical_increment: float :type azimuth: Azimuth, in decimal degrees. :type azimuth: float :returns: The point at the given distances. :rtype: Instance of :class:`Point` """ lon, lat = geodetic.point_at(self.longitude, self.latitude, azimuth, horizontal_distance) return Point(lon, lat, self.depth + vertical_increment)
[docs] def azimuth(self, point): """ Compute the azimuth (in decimal degrees) between this point and the given point. :param point: Destination point. :type point: Instance of :class:`Point` :returns: The azimuth, value in a range ``[0, 360)``. :rtype: float """ return geodetic.azimuth(self.longitude, self.latitude, point.longitude, point.latitude)
[docs] def distance(self, point): """ Compute the distance (in km) between this point and the given point. Distance is calculated using pythagoras theorem, where the hypotenuse is the distance and the other two sides are the horizontal distance (great circle distance) and vertical distance (depth difference between the two locations). :param point: Destination point. :type point: Instance of :class:`Point` :returns: The distance. :rtype: float """ return geodetic.distance(self.longitude, self.latitude, self.depth, point.longitude, point.latitude, point.depth)
[docs] def distance_to_mesh(self, mesh, with_depths=True): """ Compute distance (in km) between this point and each point of ``mesh``. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to calculate distance to. :param with_depths: If ``True`` (by default), distance is calculated between actual point and the mesh, geodetic distance of projections is combined with vertical distance (difference of depths). If this is set to ``False``, only geodetic distance between projections is calculated. :returns: Numpy array of floats of the same shape as ``mesh`` with distance values in km in respective indices. """ if with_depths: if mesh.depths is None: mesh_depths = numpy.zeros_like(mesh.lons) else: mesh_depths = mesh.depths return geodetic.distance(self.longitude, self.latitude, self.depth, mesh.lons, mesh.lats, mesh_depths) else: return geodetic.geodetic_distance(self.longitude, self.latitude, mesh.lons, mesh.lats)
def __str__(self): """ >>> str(Point(1, 2, 3)) '<Latitude=2.000000, Longitude=1.000000, Depth=3.0000>' >>> str(Point(1.0 / 3.0, -39.999999999, 1.6666666666)) '<Latitude=-40.000000, Longitude=0.333333, Depth=1.6667>' """ return "<Latitude=%.6f, Longitude=%.6f, Depth=%.4f>" % ( self.latitude, self.longitude, self.depth ) def __repr__(self): """ >>> str(Point(1, 2, 3)) == repr(Point(1, 2, 3)) True """ return self.__str__() def __eq__(self, other): """ >>> Point(1e-4, 1e-4) == Point(0, 0) False >>> Point(1e-6, 1e-6) == Point(0, 0) True >>> Point(0, 0, 1) == Point(0, 0, 0) False >>> Point(4, 5, 1e-3) == Point(4, 5, 0) True >>> Point(-180 + 1e-7, 0) == Point(180 - 1e-7, 0) True """ if other is None: return False return abs(self.distance(other)) <= self.EQUALITY_DISTANCE def __ne__(self, other): return not self.__eq__(other)
[docs] def on_surface(self): """ Check if this point is defined on the surface (depth is 0.0). :returns bool: True if this point is on the surface, false otherwise. """ return self.depth == 0.0
[docs] def equally_spaced_points(self, point, distance): """ Compute the set of points equally spaced between this point and the given point. :param point: Destination point. :type point: Instance of :class:`Point` :param distance: Distance between points (in km). :type distance: float :returns: The list of equally spaced points. :rtype: list of :class:`Point` instances """ lons, lats, depths = geodetic.intervals_between( self.longitude, self.latitude, self.depth, point.longitude, point.latitude, point.depth, distance) return [Point(lons[i], lats[i], depths[i]) for i in range(len(lons))]
[docs] def to_polygon(self, radius): """ Create a circular polygon with specified radius centered in the point. :param radius: Required radius of a new polygon, in km. :returns: Instance of :class:`~openquake.hazardlib.geo.polygon.Polygon` that approximates a circle around the point with specified radius. """ assert radius > 0 # avoid circular imports from openquake.hazardlib.geo.polygon import Polygon # get a projection that is centered in the point proj = geo_utils.OrthographicProjection( self.longitude, self.longitude, self.latitude, self.latitude) # create a shapely object from a projected point coordinates, # which are supposedly (0, 0) point = shapely.geometry.Point(*proj(self.longitude, self.latitude)) # extend the point to a shapely polygon using buffer() # and create openquake.hazardlib.geo.polygon.Polygon object from it return Polygon._from_2d(point.buffer(radius), proj)
[docs] def closer_than(self, mesh, radius): """ Check for proximity of points in the ``mesh``. :param mesh: :class:`openquake.hazardlib.geo.mesh.Mesh` instance. :param radius: Proximity measure in km. :returns: Numpy array of boolean values in the same shape as the mesh coordinate arrays with ``True`` on indexes of points that are not further than ``radius`` km from this point. Function :func:`~openquake.hazardlib.geo.geodetic.distance` is used to calculate distances to points of the mesh. Points of the mesh that lie exactly ``radius`` km away from this point also have ``True`` in their indices. """ dists = geodetic.distance(self.longitude, self.latitude, self.depth, mesh.lons, mesh.lats, 0 if mesh.depths is None else mesh.depths) return dists <= radius
[docs] @classmethod def from_vector(cls, vector): """ Create a point object from a 3d vector in Cartesian space. :param vector: Tuple, list or numpy array of three float numbers representing point coordinates in Cartesian 3d space. :returns: A :class:`Point` object created from those coordinates. """ return cls(*geo_utils.cartesian_to_spherical(vector))