# Geographic primitives¶

Various primitives are needed to define `seismic sources` location and shape. Those implemented include `point`, `line` and `polygon`.

`Mesh` objects are used only internally.

## Point¶

Module `openquake.hazardlib.geo.point` defines `Point`.

class `openquake.hazardlib.geo.point.``Point`(longitude, latitude, depth=0.0)[source]

This class represents a geographical point in terms of longitude, latitude, and depth (with respect to the Earth surface).

Parameters: longitude (float) – Point longitude, in decimal degrees. latitude (float) – Point latitude, in decimal degrees. depth (float) – Point depth (default to 0.0), in km. Depth > 0 indicates a point below the earth surface, and depth < 0 above the earth surface.
`EQUALITY_DISTANCE` = 0.001

The distance between two points for them to be considered equal, in km.

`azimuth`(point)[source]

Compute the azimuth (in decimal degrees) between this point and the given point.

Parameters: point (Instance of `Point`) – Destination point. The azimuth, value in a range `[0, 360)`. float
`closer_than`(mesh, radius)[source]

Check for proximity of points in the `mesh`.

Parameters: mesh – `openquake.hazardlib.geo.mesh.Mesh` instance. radius – Proximity measure in km. 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 `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.
`distance`(point)[source]

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).

Parameters: point (Instance of `Point`) – Destination point. The distance. float
`distance_to_mesh`(mesh, with_depths=True)[source]

Compute distance (in km) between this point and each point of `mesh`.

Parameters: mesh – `Mesh` of points to calculate distance to. 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. Numpy array of floats of the same shape as `mesh` with distance values in km in respective indices.
`equally_spaced_points`(point, distance)[source]

Compute the set of points equally spaced between this point and the given point.

Parameters: point (Instance of `Point`) – Destination point. distance (float) – Distance between points (in km). The list of equally spaced points. list of `Point` instances
classmethod `from_vector`(vector)[source]

Create a point object from a 3d vector in Cartesian space.

Parameters: vector – Tuple, list or numpy array of three float numbers representing point coordinates in Cartesian 3d space. A `Point` object created from those coordinates.
`on_surface`()[source]

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.
`point_at`(horizontal_distance, vertical_increment, azimuth)[source]

Compute the point with given horizontal, vertical distances and azimuth from this point.

Parameters: horizontal_distance (float) – Horizontal distance, in km. vertical_increment (float) – Vertical increment, in km. When positive, the new point has a greater depth. When negative, the new point has a smaller depth. The point at the given distances. Instance of `Point`
`to_polygon`(radius)[source]

Create a circular polygon with specified radius centered in the point.

Parameters: radius – Required radius of a new polygon, in km. Instance of `Polygon` that approximates a circle around the point with specified radius.
`wkt2d`

Generate WKT (Well-Known Text) to represent this point in 2 dimensions (ignoring depth).

`x`

Alias for .longitude

`y`

Alias for .latitude

`z`

Alias for .depth

## Line¶

Module `openquake.hazardlib.geo.line` defines `Line`.

class `openquake.hazardlib.geo.line.``Line`(points)[source]

This class represents a geographical line, which is basically a sequence of geographical points.

A line is defined by at least one point.

Parameters: points (list of `Point` instances) – The sequence of points defining this line.
`average_azimuth`()[source]

Calculate and return weighted average azimuth of all line’s segments in decimal degrees.

Uses formula from http://en.wikipedia.org/wiki/Mean_of_circular_quantities

```>>> from openquake.hazardlib.geo.point import Point as P
>>> str(Line([P(0, 0), P(1e-5, 1e-5)]).average_azimuth())
'45.0'
>>> str(Line([P(0, 0), P(0, 1e-5), P(1e-5, 1e-5)]).average_azimuth())
'45.0'
>>> line = Line([P(0, 0), P(-2e-5, 0), P(-2e-5, 1.154e-5)])
>>> '%.1f' % line.average_azimuth()
'300.0'
```
`get_length`()[source]

Calculate and return the length of the line as a sum of lengths of all its segments.

Returns: Total length in km.
`on_surface`()[source]

Check if this line is defined on the surface (i.e. all points are on the surfance, depth=0.0).

Returns bool: True if this line is on the surface, false otherwise.
`resample`(section_length)[source]

Resample this line into sections.

The first point in the resampled line corresponds to the first point in the original line.

Starting from the first point in the original line, a line segment is defined as the line connecting the last point in the resampled line and the next point in the original line. The line segment is then split into sections of length equal to `section_length`. The resampled line is obtained by concatenating all sections.

The number of sections in a line segment is calculated as follows: `round(segment_length / section_length)`.

Note that the resulting line has a length that is an exact multiple of `section_length`, therefore its length is in general smaller or greater (depending on the rounding) than the length of the original line.

For a straight line, the difference between the resulting length and the original length is at maximum half of the `section_length`. For a curved line, the difference my be larger, because of corners getting cut.

Parameters: section_length (float) – The length of the section, in km. A new line resampled into sections based on the given length. An instance of `Line`
`resample_to_num_points`(num_points)[source]

Resample the line to a specified number of points.

Parameters: num_points – Integer number of points the resulting line should have. A new line with that many points as requested.

## Polygon¶

class `openquake.hazardlib.geo.polygon.``Polygon`(points)[source]

Polygon objects represent an area on the Earth surface.

Parameters: points – The list of `Point` objects defining the polygon vertices. The points are connected by great circle arcs in order of appearance. Polygon segment should not cross another polygon segment. At least three points must be defined. ValueError – If `points` contains less than three unique points or if polygon perimeter intersects itself.
`dilate`(dilation)[source]

Extend the polygon to a specified buffer distance.

Note

In extreme cases where dilation of a polygon creates holes, thus resulting in a multi-polygon, we discard the holes and simply return the ‘exterior’ of the shape.

Parameters: dilation – Distance in km to extend polygon borders to. New `Polygon` object with (in general) more vertices and border that is approximately `dilation` km far (measured perpendicularly to edges and circularly to vertices) from the border of original polygon.
`discretize`(mesh_spacing)[source]

Get a mesh of uniformly spaced points inside the polygon area with distance of `mesh_spacing` km between.

Returns: An instance of `Mesh` that holds the points data. Mesh is created with no depth information (all the points are on the Earth surface).
`intersects`(mesh)[source]

Check for intersection with each point of the `mesh`.

Mesh coordinate values are in decimal degrees.

Parameters: mesh – `openquake.hazardlib.geo.mesh.Mesh` instance. Numpy array of bool values in the same shapes in the input coordinate arrays with `True` on indexes of points that lie inside the polygon or on one of its edges and `False` for points that neither lie inside nor touch the boundary.
`wkt`

Generate WKT (Well-Known Text) to represent this polygon.

`openquake.hazardlib.geo.polygon.``UPSAMPLING_STEP_KM` = 100

Polygon upsampling step for long edges, in kilometers. See `get_resampled_coordinates()`.

`openquake.hazardlib.geo.polygon.``get_resampled_coordinates`(lons, lats)[source]

Resample polygon line segments and return the coordinates of the new vertices. This limits distortions when projecting a polygon onto a spherical surface.

Parameters define longitudes and latitudes of a point collection in the form of lists or numpy arrays.

Returns: A tuple of two numpy arrays: longitudes and latitudes of resampled vertices.

## Meshes¶

Module `openquake.hazardlib.geo.mesh` defines classes `Mesh` and its subclass `RectangularMesh`.

### Simple mesh¶

class `openquake.hazardlib.geo.mesh.``Mesh`(lons, lats, depths=None)[source]

Mesh object represent a collection of points and provides the most efficient way of keeping those collections in memory.

Parameters: lons – A numpy array of longitude values of points. Array may be of arbitrary shape. lats – Numpy array of latitude values. The array must be of the same shape as `lons`. depths – Either `None`, which means that all points the mesh consists of are lying on the earth surface (have zero depth) or numpy array of the same shape as previous two.

Mesh object can also be created from a collection of points, see `from_points_list()`.

`DIST_TOLERANCE` = 0.005

Tolerance level to be used in various spatial operations when approximation is required – set to 5 meters.

classmethod `from_points_list`(points)[source]

Create a mesh object from a collection of points.

Parameters: point – List of `Point` objects. An instance of `Mesh` with one-dimensional arrays of coordinates from `points`.
`get_closest_points`(mesh)[source]

Find closest point of this mesh for each one in `mesh`.

Returns: `Mesh` object of the same shape as `mesh` with closest points from this one at respective indices.

This method is in general very similar to `get_min_distance()` and uses the same `openquake.hazardlib.geo.geodetic.min_distance()` internally.

`get_convex_hull`()[source]

Get a convex polygon object that contains projections of all the points of the mesh.

Returns: Instance of `openquake.hazardlib.geo.polygon.Polygon` that is a convex hull around all the points in this mesh. If the original mesh had only one point, the resulting polygon has a square shape with a side length of 10 meters. If there were only two points, resulting polygon is a stripe 10 meters wide.
`get_distance_matrix`()[source]

Compute and return distances between each pairs of points in the mesh.

This method requires that all the points lie on Earth surface (have zero depth) and coordinate arrays are one-dimensional.

Warning

Because of its quadratic space and time complexity this method is safe to use for meshes of up to several thousand points. For mesh of 10k points it needs ~800 Mb for just the resulting matrix and four times that much for intermediate storage.

Returns: Two-dimensional numpy array, square matrix of distances. The matrix has zeros on main diagonal and positive distances in kilometers on all other cells. That is, value in cell (3, 5) is the distance between mesh’s points 3 and 5 in km, and it is equal to value in cell (5, 3).
`get_joyner_boore_distance`(mesh)[source]

Compute and return Joyner-Boore distance to each point of `mesh`. Point’s depth is ignored.

See `openquake.hazardlib.geo.surface.base.BaseQuadrilateralSurface.get_joyner_boore_distance()` for definition of this distance.

Returns: numpy array of distances in km of the same shape as `mesh`. Distance value is considered to be zero if a point lies inside the polygon enveloping the projection of the mesh or on one of its edges.
`get_min_distance`(mesh)[source]

Compute and return the minimum distance from the mesh to each point in another mesh.

Returns: numpy array of distances in km of the same shape as `mesh`.

Method doesn’t make any assumptions on arrangement of the points in either mesh and instead calculates the distance from each point of this mesh to each point of the target mesh and returns the lowest found for each.

`shape`

Return the shape of this mesh.

Returns tuple: The shape of this mesh as (rows, columns)

### Rectangular mesh¶

class `openquake.hazardlib.geo.mesh.``RectangularMesh`(lons, lats, depths=None)[source]

A specification of `Mesh` that requires coordinate numpy-arrays to be two-dimensional.

Rectangular mesh is meant to represent not just an unordered collection of points but rather a sort of table of points, where index of the point in a mesh is related to it’s position with respect to neighbouring points.

classmethod `from_points_list`(points)[source]

Create a rectangular mesh object from a list of lists of points. Lists in a list are supposed to have the same length.

Parameters: point – List of lists of `Point` objects.
`get_cell_dimensions`()[source]

Calculate centroid, width, length and area of each mesh cell.

Returns: Tuple of four elements, each being 2d numpy array. Each array has both dimensions less by one the dimensions of the mesh, since they represent cells, not vertices. Arrays contain the following cell information: centroids, 3d vectors in a Cartesian space, length (size along row of points) in km, width (size along column of points) in km, area in square km.
`get_mean_inclination_and_azimuth`()[source]

Calculate weighted average inclination and azimuth of the mesh surface.

Returns: Tuple of two float numbers: inclination angle in a range [0, 90] and azimuth in range [0, 360) (in decimal degrees).

The mesh is triangulated, the inclination and azimuth for each triangle is computed and average values weighted on each triangle’s area are calculated. Azimuth is always defined in a way that inclination angle doesn’t exceed 90 degree.

`get_mean_width`()[source]

Calculate and return (weighted) mean width (km) of a mesh surface.

The length of each mesh column is computed (summing up the cell widths in a same column), and the mean value (weighted by the mean cell length in each column) is returned.

`get_middle_point`()[source]

Return the middle point of the mesh.

Returns: An instance of `Point`.

The middle point is taken from the middle row and a middle column of the mesh if there are odd number of both. Otherwise the geometric mean point of two or four middle points.

`triangulate`()[source]

Convert mesh points to vectors in Cartesian space.

Returns: Tuple of four elements, each being 2d numpy array of 3d vectors (the same structure and shape as the mesh itself). Those arrays are: points vectors, vectors directed from each point (excluding the last column) to the next one in a same row →, vectors directed from each point (excluding the first row) to the previous one in a same column ↑, vectors pointing from a bottom left point of each mesh cell to top right one ↗. So the last three arrays of vectors allow to construct triangles covering the whole mesh.