Source code for mpas_tools.transects

import numpy
from shapely.geometry import LineString, Point


[docs]def subdivide_great_circle(x, y, z, maxRes, earthRadius): """ Subdivide each segment of the transect so the horizontal resolution approximately matches the requested resolution Uses a formula for interpolating unit vectors on the sphere from https://en.wikipedia.org/wiki/Slerp Parameters ---------- x : numpy.ndarray The Cartesian x coordinate of a transect, where the number of segments is ``len(x) - 1``. ``x``, ``y`` and ``z`` are of the same length. y : numpy.ndarray The Cartesian y coordinate of the transect z : numpy.ndarray The Cartesian z coordinate of the transect maxRes : float The maximum allowed spacing in m after subdivision earthRadius : float The radius of the Earth in m Returns ------- xOut : numpy.ndarray The Cartesian x values of the transect subdivided into segments with segment length at most ``maxRes``. All the points in ``x``, ``y`` and ``z`` are guaranteed to be included. yOut : numpy.ndarray The Cartesian y values of the subdivided transect zOut : numpy.ndarray The Cartesian y values of the subdivided transect dIn : numpy.ndarray The distance along the transect before subdivision dOut : numpy.ndarray The distance along the transect after subdivision """ angularDistance = angular_distance(x=x, y=y, z=z) dx = angularDistance * earthRadius nSegments = numpy.maximum( (dx / maxRes + 0.5).astype(int), 1) dIn = numpy.zeros(x.shape) dIn[1:] = numpy.cumsum(dx) frac = [] outIndices = [] delta = [] for index in range(len(dIn) - 1): n = nSegments[index] frac.extend(numpy.arange(0, n)/n) outIndices.extend(index*numpy.ones(n, int)) delta.extend(angularDistance[index]*numpy.ones(n)) frac.append(1.) outIndices.append(len(dIn) - 2) delta.append(angularDistance[-1]) frac = numpy.array(frac) delta = numpy.array(delta) outIndices = numpy.array(outIndices) a = numpy.ones(delta.shape) b = numpy.zeros(delta.shape) mask = delta > 0. denom = 1./numpy.sin(delta[mask]) a[mask] = denom*numpy.sin((1.-frac[mask])*delta[mask]) b[mask] = denom*numpy.sin(frac[mask]*delta[mask]) xOut = a*x[outIndices] + b*x[outIndices+1] yOut = a*y[outIndices] + b*y[outIndices+1] zOut = a*z[outIndices] + b*z[outIndices+1] dOut = (-frac + 1.)*dIn[outIndices] + frac*dIn[outIndices+1] return xOut, yOut, zOut, dIn, dOut
[docs]def cartesian_to_great_circle_distance(x, y, z, earth_radius): """ Cartesian transect points to great-circle distance Parameters ---------- x : numpy.ndarray The Cartesian x coordinate of a transect y : numpy.ndarray The Cartesian y coordinate of the transect z : numpy.ndarray The Cartesian z coordinate of the transect earth_radius : float The radius of the earth Returns ------- distance : numpy.ndarray The distance along the transect """ distance = numpy.zeros(x.shape) for segIndex in range(len(x)-1): transectv0 = Vector(x[segIndex], y[segIndex], z[segIndex]) transectv1 = Vector(x[segIndex+1], y[segIndex+1], z[segIndex+1]) distance[segIndex+1] = distance[segIndex] + \ earth_radius*angular_distance(first=transectv0, second=transectv1) return distance
[docs]def subdivide_planar(x, y, maxRes): """ Subdivide each segment of the transect so the horizontal resolution approximately matches the requested resolution Uses a formula for interpolating unit vectors on the sphere from https://en.wikipedia.org/wiki/Slerp Parameters ---------- x : numpy.ndarray The planar x coordinate of a transect, where the number of segments is ``len(x) - 1`` y : numpy.ndarray The planar y coordinates of the transect, the same length as ``x`` maxRes : float The maximum allowed spacing in m after subdivision Returns ------- xOut : numpy.ndarray The x coordinate of the transect, subdivided into segments with length at most ``maxRes``. All the points in ``x`` are guaranteed to be included. yOut : numpy.ndarray The y coordinate of the transect. All the points in ``y`` are guaranteed to be included. dIn : numpy.ndarray The distance along the transect before subdivision dOut : numpy.ndarray The distance along the transect after subdivision """ dx = numpy.zeros(len(x)-1) for index in range(len(x)-1): start = Point(x[index], y[index]) end = Point(x[index+1], y[index+1]) segment = LineString([start, end]) dx[index] = segment.length nSegments = numpy.maximum( (dx / maxRes + 0.5).astype(int), 1) dIn = numpy.zeros(x.shape) dIn[1:] = numpy.cumsum(dx) frac = [] outIndices = [] for index in range(len(dIn) - 1): n = nSegments[index] frac.extend(numpy.arange(0, n)/n) outIndices.extend(index*numpy.ones(n, int)) frac.append(1.) outIndices.append(len(dIn) - 2) frac = numpy.array(frac) outIndices = numpy.array(outIndices) xOut = (-frac + 1.)*x[outIndices] + frac*x[outIndices+1] yOut = (-frac + 1.)*y[outIndices] + frac*y[outIndices+1] dOut = (-frac + 1.)*dIn[outIndices] + frac*dIn[outIndices+1] return xOut, yOut, dIn, dOut
[docs]def lon_lat_to_cartesian(lon, lat, earth_radius, degrees): """Convert from lon/lat to Cartesian x, y, z""" if degrees: lon = numpy.deg2rad(lon) lat = numpy.deg2rad(lat) x = earth_radius * numpy.cos(lat) * numpy.cos(lon) y = earth_radius * numpy.cos(lat) * numpy.sin(lon) z = earth_radius * numpy.sin(lat) return x, y, z
[docs]def cartesian_to_lon_lat(x, y, z, earth_radius, degrees): """Convert from Cartesian x, y, z to lon/lat""" lon = numpy.arctan2(y, x) lat = numpy.arcsin(z/earth_radius) if degrees: lon = numpy.rad2deg(lon) lat = numpy.rad2deg(lat) return lon, lat
[docs]def angular_distance(x=None, y=None, z=None, first=None, second=None): """ Compute angular distance between points on the sphere, following: https://en.wikipedia.org/wiki/Great-circle_distance Parameters ---------- x : numpy.ndarray, optional The Cartesian x coordinate of a transect, where the number of segments is ``len(x) - 1``. ``x``, ``y`` and ``z`` are of the same length and all must be present if ``first`` and ``second`` are not provided. y : numpy.ndarray, optional The Cartesian y coordinate of the transect z : numpy.ndarray, optional The Cartesian z coordinate of the transect first : mpas_tools.transect.Vector, optional The start points of each segment of the transect, where the ``x``, ``y``, and ``z`` attributes of each vector are numpy.ndarray objects. second : mpas_tools.transect.Vector, optional The end points of each segment of the transect Returns ------- angularDistance : numpy.ndarray The angular distance (in radians) between segments of the transect. """ if first is None or second is None: first = Vector(x[0:-1], y[0:-1], z[0:-1]) second = Vector(x[1:], y[1:], z[1:]) angularDistance = numpy.arctan2(_mag(_cross(first, second)), _dot(first, second)) return angularDistance
[docs]def intersects(a1, a2, b1, b2): """ Based on https://stackoverflow.com/a/26669130/7728169 Determine if the great circle arc from ``a1`` to ``a2`` intersects that from ``b1`` to ``b2``. Parameters ---------- a1 : mpas_tools.transects.Vector Cartesian coordinates of the end point of a great circle arc. The types of the attributes ``x``, ``y``, and ``z`` must either be ``numpy.arrays`` of identical size for all 4 vectors (in which case intersections are found element-wise), or scalars for at least one of either ``a1`` and ``a2`` or ``b1`` and ``b2``. a2 : mpas_tools.transects.Vector Cartesian coordinates of the other end point of a great circle arc. b1 : mpas_tools.transects.Vector Cartesian coordinates of an end point of a second great circle arc. b2 : mpas_tools.transects.Vector Cartesian coordinates of the other end point of the second great circle arc. Returns ------- intersect : numpy.ndarray A boolean array of the same size as ``a1`` and ``a2`` or ``b1`` and ``b2``, whichever is greater, indicating if the particular pair of arcs intersects """ return numpy.logical_and(_straddles(a1, a2, b1, b2), _straddles(b1, b2, a1, a2))
[docs]def intersection(a1, a2, b1, b2): """ Based on https://stackoverflow.com/a/26669130/7728169 Find the intersection point as a unit vector between great circle arc from ``a1`` to ``a2`` and from ``b1`` to ``b2``. The arcs should have already have been found to intersect by calling ``intersects()`` Parameters ---------- a1 : mpas_tools.transects.Vector Cartesian coordinates of the end point of a great circle arc. The types of the attributes ``x``, ``y``, and ``z`` must either be ``numpy.arrays`` of identical size for all 4 vectors (in which case intersections are found element-wise), or scalars for at least one of either ``a1`` and ``a2`` or ``b1`` and ``b2``. a2 : mpas_tools.transects.Vector Cartesian coordinates of the other end point of a great circle arc. b1 : mpas_tools.transects.Vector Cartesian coordinates of an end point of a second great circle arc. b2 : mpas_tools.transects.Vector Cartesian coordinates of the other end point of the second great circle arc. Returns ------- points : mpas_tools.transects.Vector An array of Cartesian points *on the unit sphere* indicating where the arcs intersect """ points = _cross(_cross(a1, a2), _cross(b1, b2)) s = numpy.sign(_det(a1, b1, b2))/_mag(points) points = Vector(s*points.x, s*points.y, s*points.z) return points
[docs]class Vector: """ A class for representing Cartesian vectors with ``x``, ``y`` and ``z`` components that are either ``float`` or ``numpy.array`` objects of identical size. """
[docs] def __init__(self, x, y, z): self.x = x self.y = y self.z = z
def _straddles(a1, a2, b1, b2): """ Based on https://stackoverflow.com/a/26669130/7728169 Determines if the great circle segment determined by (a1, a2) straddles the great circle determined by (b1, b2) Parameters ---------- a1, a2, b1, b2 : Vector Cartesian coordinates of the end points of two great circle arcs. The types of the attributes ``x``, ``y``, and ``z`` must either be ``numpy.arrays`` of identical size for all 4 vectors (in which case intersections are found element-wise), or scalars for at least one of either the ``a``s or the ``b``s. Returns ------- straddle : numpy.ndarray A boolean array of the same size as the ``a``s or the ``b``s, whichever is greater, indicating if the great circle segment determined by (a1, a2) straddles the great circle determined by (b1, b2) """ return _det(a1, b1, b2) * _det(a2, b1, b2) < 0 def _dot(v1, v2): """The dot product between two ``Vector`` objects ``v1`` and ``v2``""" return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z def _cross(v1, v2): """The cross product between two ``Vector`` objects ``v1`` and ``v2``""" return Vector(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x) def _det(v1, v2, v3): """The determinant of the matrix defined by the three ``Vector`` objects""" return _dot(v1, _cross(v2, v3)) def _mag(v): """The magnitude of the ``Vector`` object ``v``""" return numpy.sqrt(_dot(v, v))