Source code for irfpy.util.tessellatedsphere

''' Tessellating a sphere surface by triangles.

For mapping data to a sphere, tessellation is a powerful way.
Usually, people work with "cartesian", the "longitude-latitude" system.
It is widely used as the binning is intuitive and straightforward (indeed not as so as it looks),
and easy to implement.  A sphere is split into, for example, 360 longitude bins
and 180 latitude bin.
However, the longitude-latitude system has weakness close to the polar regions.
The size of the bin depends on cosine of the longitude, which potentially produces a lot of problems.

Here trianglar tessellation comes into the game.

Here I implemented

- :class:`RecurrentTessellatedUnitSphere`: An implementation of the tessellated sphere, a spphere represented by
  multiple triangles with almost the same area size.
- :class:`TessellatedSphere`: A simple implementation of the tessellated sphere, a sphere represented by
  multiple triangles with almost the same area size.  However, this implementation is for reference and legacy,
  so that use the :class:`RecurrentTessellatedUnitSphere`

What we can do

- Obtain multiple triangles (three vortex coordinates) that forms a sphere
- Obtain which triangle contain a given vector (conversion from a coordinates to a triangle)
- (Advanced) Relation between triangles (links, i.e., which triangles are the neighbor of the triangle of interest?)

.. codeauthor:: Yoshifumi Futaana
'''
import time
import logging

import numpy as np
from irfpy.util import triangle
from irfpy.util import icosahedron

from irfpy.util import vector3d


[docs]class TessellatedSphere: ''' A tessellated sphere. This class represents a tessellated sphere. .. warning:: This class is for reference implementation using legacy irfpy.util functionality. Use :class:`RecurrentTessellatedUnitSphere` for practical use. Only a single parameter *ndim* is taken; here *ndim* is the degree of resolution. - If *ndim* is 0, it is just a regular icosahedron. - If *ndim* is 1, each triangle of the icosahedron (ndim=0) is split into 4 by choosing the central point of the triangle edges mapped to the sphere surface. This means that the number of the triangles becomes 80. - For *ndim* of N, each triangle forming *ndim* level (N-1) of tessellated sphere is split into 4. This results in the N-tessellated sphere has 4 times higher resolution (number of triangles) than the (N-1)-tessellated sphere. ''' def __init__(self, ndim=0, radius=1.0): ''' The constractor. :keyword ndim: The degree of tessellation of the sphere. Default is 0, the icosahedron. :type ndim: ``int`` :keyword radius: Radius of the sphere. Default is 1. :type radius: ``float`` You can create the tessellated sphere with specifying the ndim, the degree of details. Here for example, the 0-level tessellated sphere has 20 triangle, the normal icoshedron. >>> ts = TessellatedSphere(ndim=0) >>> tris = ts.getTriangles() # You can get triangles. >>> print(len(tris)) 20 >>> print(tris[0]) Triangle(Vector3d( 0, 0, 1 ) Vector3d( 0.894427, 0, 0.447214 ) Vector3d( 0.276393, 0.850651, 0.447214 )) :meth:`getTriangles` gives back all the triangles. Triangle is implemented as :class:`irfpy.util.triangle.Triangle` but the triangle class has been deprecated. Further detailed tessellated spheres can be obtained as follows. >>> ts = TessellatedSphere(ndim=1) >>> tris = ts.getTriangles() >>> print(len(tris)) 80 >>> ts = TessellatedSphere(ndim=2) >>> tris = ts.getTriangles() >>> print(len(tris)) 320 >>> ts = TessellatedSphere(ndim=3) >>> tris = ts.getTriangles() >>> print(len(tris)) 1280 The 3rd-level tessellated sphere has 1280 triangles. The total area is almost the same as 4-pi. >>> areas = [_tri.getArea() for _tri in tris] >>> print('{:.3f}'.format(np.sum(areas))) 12.506 .. note:: **Note for developer** Internally, the radius of the circumscribe sphere is 1.0, i.e. the unit sphere. This means that the radius is only the parameter when users get :class:`Triangle`. If you want a tessellated sphere with radius of 7, you can specify the radius in the constructor. >>> ts3 = TessellatedSphere(ndim=3, radius=7) ''' if ndim is None: self.ndim = 0 else: self.ndim = ndim self.radius = radius self.tris = self._makeTessellation(self.ndim) def _makeTessellation(self, ndim): ''' Calculate the tessellation. (Called internally only) :param ndim: Tessellation level. :type ndim: int :returns: A list of tessellated triangle. :rtype: List of Triangle ''' # Initialize ic = icosahedron.RegularIcosahedron() tris = ic.getTriangles() # tris is a tuple of Triangle instances. for i in range(ndim): tris = self._separateTriangleList(tris) return tris def _separateTriangleList(self, tris): ''' List of triangles are splitted. Each the trinagle in the given list is split into 4, and the resulting triangles are returned as a list. Only called internally. ''' newtris = [] for tri in tris: ntris = self._separateTriangle(tri) for ntri in ntris: newtris.append(ntri) return newtris @classmethod def _separateTriangle(self, tri, strictness=1e-7): r''' A triangle is split into four triangles. :Parameters: tri : Triangle A triangle to be separated. Given triangles should have been normalized for this method. Check will be done in the beginning of the method. *strictness* parameter can change the behaviour. strictness : float If any of the vertices has lengthSquared > 1 +- *strictness*, RuntimeError is raised. :Returns: t0, t1, t2, t3 : Triangle Split triangles. Splitting the triangle is like this :: v0 ---- v5 ---- v2 \ / \ / v3 --- v4 \ / v1 The triangle ``(v0, v1, v2)`` is split into four triangles of - ``(v0, v3, v5)`` - ``(v1, v4, v3)`` - ``(v2, v5, v4)`` - ``(v3, v4, v5)`` Here, v3, v4, and v5 is the center point of v0, v1, and v2, and then mapped to the surface. >>> v0 = vector3d.Vector3d(1, 5, -2) >>> v0.normalize() >>> v1 = vector3d.Vector3d(3, 2, -5) >>> v1.normalize() >>> v2 = vector3d.Vector3d(2, -1, 3) >>> v2.normalize() >>> t = triangle.Triangle(v0, v1, v2) >>> t0, t1, t2, t3 = TessellatedSphere._separateTriangle(t) >>> v00, v01, v02 = t0.getVertices() >>> v00.sub(v0) >>> v00.lengthSquared() < 1e-7 True Called only internally. ''' v0, v1, v2 = tri.getVertices() if abs(v0.lengthSquared() - 1) > strictness: raise RuntimeError('Vector v0 looks non-normalized. (Leng^2 = %g)' % v0.lengthSquared()) if abs(v1.lengthSquared() - 1) > strictness: raise RuntimeError('Vector v1 looks non-normalized. (Leng^2 = %g)' % v1.lengthSquared()) if abs(v2.lengthSquared() - 1) > strictness: raise RuntimeError('Vector v2 looks non-normalized. (Leng^2 = %g)' % v2.lengthSquared()) v3 = vector3d.Vector3d(v0.x + v1.x, v0.y + v1.y, v0.z + v1.z) v3.normalize() v4 = vector3d.Vector3d(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z) v4.normalize() v5 = vector3d.Vector3d(v2.x + v0.x, v2.y + v0.y, v2.z + v0.z) v5.normalize() t0 = triangle.Triangle(v0, v3, v5) t1 = triangle.Triangle(v1, v4, v3) t2 = triangle.Triangle(v2, v5, v4) t3 = triangle.Triangle(v3, v4, v5) return (t0, t1, t2, t3)
[docs] def getTriangles(self): ''' Returns the tuple of the triangles that tessellate the sphere. ''' tris = [] for tri in self.tris: v0, v1, v2 = tri.getVertices() v0.scale(self.radius) v1.scale(self.radius) v2.scale(self.radius) tris.append(triangle.Triangle(v0, v1, v2)) return tuple(tris)
[docs]class RecurrentTessellatedUnitSphere: r''' A tessellated sphere, a list of triangles that forms the surface of the sphere. This class represents a tessellated sphere. Only a single parameter *ndim* is taken; here *ndim* is the degree of resolution. - If *ndim* is 0, it is just a regular icosahedron. - If *ndim* is 1, each triangle of the icosahedron (ndim=0) is split into 4 by choosing the central point of the triangle edges mapped to the sphere surface. This means that the number of the triangles becomes 80. - For *ndim* of N, each triangle forming *ndim* level (N-1) of tessellated sphere is split into 4. This results in the N-tessellated sphere has 4 times higher resolution (number of triangles) than the (N-1)-tessellated sphere. The number of the triangles are :math:`80\times 4^N`. - ndim=3 gives 1280, which may be good for testing. - ndim=5 gives 20480, which may be good for normal production - ndim=6 gives 81920, which might give the good final resolution that can be handled for usual scientific purpose. Note that 180x360 (1 degree grid) for classical use is 67800 bins. - ndim=7 gives 327680 bins. This may work, but each processing may time a while (seconds to minutes). Some functions are very slow. This indicates that this is almost the hightest resolution that the current implementation could handle techincally. ''' logger = logging.getLogger(__name__ + '.RecurrentTessellatedUnitSphere') def __init__(self, ndim=0): ''' The constractor. :keyword ndim: The degree of tessellation of the sphere. Default is 0, the icosahedron. :type ndim: ``int`` You can create the tessellated sphere with specifying the ndim, the degree of details. Here for example, the 0-level tessellated sphere has 20 triangle, the normal icoshedron. >>> ts = RecurrentTessellatedUnitSphere(ndim=0) >>> tris = ts.getTriangles() # You can get triangles. >>> print(len(tris)) 20 Further detailed tessellated spheres can be obtained as follows. >>> ts = RecurrentTessellatedUnitSphere(ndim=1) >>> tris = ts.getTriangles() >>> print(len(tris)) 80 >>> print(tris[0]) NpTriangle([0. 0. 1.] [0.52573111 0. 0.85065081] [0.16245985 0.5 0.85065081]) :meth:`getTriangles` gives back all the triangles. Triangle is implemented as :class:`irfpy.util.triangle.NpTriangle`. You may want to get the full numpy array for triangles. >>> triangle_vertexes = ts.getTriangleVertexes() >>> print(triangle_vertexes.shape) (80, 3, 3) >>> print(triangle_vertexes[0, 0]) # Triangle #0, Vertex #0 [0. 0. 1.] >>> print(triangle_vertexes[0, 1]) # Triangle #0, Vertex #1 [0.52573111 0. 0.85065081] >>> print(triangle_vertexes[0, 2]) # Triangle #0, Vertex #2 [0.16245985 0.5 0.85065081] Instead, you also can get the longitude and latitude for all the triangles. >>> lons, lats = ts.getTriangleLonLat() >>> print(lons.shape) (80, 3) >>> print(lons[30, 1]) # Longitude of the triangle #30, the vertex 1. 108.0 You may find the triangle corresponding to the given position. For example, let's see the longitude, latitude as >>> lon = 20 >>> lat = 67.8 Then, you can find the corresponding triangle >>> ts = RecurrentTessellatedUnitSphere(ndim=5) >>> lonlist, latlist = ts.getTriangleLonLat() >>> tindex = ts.indices_lonlat(lon, lat) >>> print(tindex) 119 >>> print(lonlist[tindex], latlist[tindex]) [16.25128128 14.93628104 20.62647783] [69.11246462 67.1720864 67.79965681] ''' t0 = time.time() self.logger.debug('Initializer: ndim={ndim}'.format(ndim=ndim)) self.ndim = ndim self.tess_tris = None ''' Tessellated triangles (ndarray of :class:`triangle.NpTriangle`) Numpy object array of :class:`triangle.NpTriangle`. (N,) shape, where N is the number of triangles. :meth:`getTriangles` is the better way to get access to this variable. ''' self.index_map = None ''' This variable describes the relation of the triangle to a triangle of the parent (ndim-1) tesseration. Since there is hyerarchic relation between the parent (ndim-1) tessellating triangle and the present (ndim) triangle (as the present triangles are the result of the split of the parent triangle), the variable keeps track of these relations. This variable is used mainly to search for the tessellated triangle corresponding to the given point, such as :meth:`indecies` or :meth:`indecies_lonlat`. The shape of this variable is (:math:`N_p`, the number of triangles of parent tessellation). If :attr:`ndim` is 0, ``None`` is set. For example, :attr:`index_map` [108] is a ndarray of indecies of triangles (for this tessellation, ndim=ndim) that are derived from the parent triangle indexed by 108. ''' self.linklist = None ''' Contains information on triangle links. :attr:`linklist` [127] will show a ndarray of indecies that shares an edge to triangle indexed by 127. If ``None``, the linklist is not initialized, so that you may need to call :meth:`createLinkList`. ''' if ndim < 0: raise ValueError('ndim=%d must be >=0' % ndim) if ndim == 0: self.parent = None self.tess_tris = self._new_triangles() else: self.parent = RecurrentTessellatedUnitSphere(ndim=ndim - 1) self.tess_tris, self.index_map = self._split_triangles(self.parent.tess_tris) self.tess_tris = np.array(self.tess_tris, dtype=np.dtype(object)) # numpy array of Triangle self.index_map = np.array(self.index_map) # (N, 4) array of index. self.norms = np.array([tri.normalVector() for tri in self.tess_tris]) # (N, 3) array self.cogs = np.array([tri.centerOfGravity() for tri in self.tess_tris]) # (N, 3) array self.areas = np.sqrt((self.norms ** 2).sum(1)) ### Normalize with efficient way. norms_t = self.norms.T # Transpose to have (3, N) array. Transpose is a light operation in numpy. norms_t /= self.areas # Then, areas (N,) shape can be in-place devided. self.norms = norms_t.T # Transpose again to have (N, 3) array. self.logger.debug('Triangle generated for ndim={ndim}. N={N}.'.format(ndim=ndim, N=len(self.tess_tris))) self.logger.debug('Elapsed time= {t} s'.format(t=time.time() - t0))
[docs] def getTriangles(self): """ Return the (reference to) list of tessellated triangles. """ return self.tess_tris
[docs] def getTriangleVertexes(self): """ Return the numpy array of the tessellated triangles. :returns: Numpy array of vertexes defined for triangles. The shape is (N, 3, 3), where N is the number of triangles (depending on the level), the next 3 is for the 3 vertexes, and the last 3 is for the (x, y, z) coordinates. >>> ts = RecurrentTessellatedUnitSphere(ndim=3) >>> ver = ts.getTriangleVertexes() >>> print(ver.shape) (1280, 3, 3) >>> print(ver[324, 1]) [0.91772872 0.37174803 0.13991924] The above example indicates that the level-3 tessellated sphere's triangle #324 has a vertex at (0.9177, 0.3717, 0.1399). """ vert = np.array([(_tris[0], _tris[1], _tris[2]) for _tris in self.tess_tris]) return vert
[docs] def getTriangleLonLat(self): """ Return the tessellating triangle longitude-latitude. Similar to :func:`getTriangleVertexes`, the function returns the triangles' longitude and latitude in degrees. :returns: Two array, (longitutde, latitude). Each array has dimension of (N, 3), where N is the number of triangles (depending on the ndim) and 3 is for the three triangle vertexes. >>> ts = RecurrentTessellatedUnitSphere(ndim=2) >>> lon, lat = ts.getTriangleLonLat() >>> print(lon.shape) (320, 3) >>> print(lat.shape) (320, 3) >>> print('{:.2f}'.format(lon[183, 2]), '{:.2f}'.format(lat[183, 2])) -117.51 -13.44 """ vert = self.getTriangleVertexes() lon = np.arctan2(vert[..., 1], vert[..., 0]) lat = np.arcsin(vert[..., 2]) return np.rad2deg(lon), np.rad2deg(lat)
def _search_index1(self, point): ''' Search the index including point. Comprehensive search. Check which center of gravity is closest. Called when ndim=0. ''' cogs = self.cogs ### Length of the cogs for each is assumed to be the same. inners = np.array([np.dot(point, cog) for cog in cogs]) self.logger.debug(inners) maxindex = np.argmax(inners) self.logger.debug('Max index = {maxindex} with value = {value}'.format( maxindex=maxindex, value=inners[maxindex])) return maxindex def _search_indecies(self, points): ''' Comprehensive search. See :meth:`indecies` for argument. (This method can work with >=3 dim array) ''' cogs = self.cogs # (20, 3) for ndim=0 inners = np.dot(points, cogs.T) # (N, ..., 3) dot (3, 20) => (N, ..., 20) maxindecies = np.argmax(inners, axis=points.ndim-1) return maxindecies
[docs] def indecies(self, points): ''' Return the indecies of the triangle IDs that contain the given points. :param points: (3,) or (n,3) shape numpy array that gives the coordinates of point(s). :returns: Scalar or n-shape numpy array. The value(s) are the indecies of triangles that contain the given points. For a single point sample >>> lo = np.deg2rad(18.7) >>> la = np.deg2rad(-24.9) >>> point = np.array([np.cos(lo) * np.cos(la), np.sin(lo) * np.cos(la), np.sin(la)]) >>> t0 = RecurrentTessellatedUnitSphere(ndim=0) >>> print(t0.indecies(point)) 14 >>> t1 = RecurrentTessellatedUnitSphere(ndim=1) >>> print(t1.indecies(point)) 58 >>> t2 = RecurrentTessellatedUnitSphere(ndim=2) >>> print(t2.indecies(point)) 235 >>> lo = np.deg2rad(-78.1) >>> la = np.deg2rad(58.2) >>> point2 = np.array([np.cos(lo) * np.cos(la), np.sin(lo) * np.cos(la), np.sin(la)]) >>> points = np.array([point, point2]) >>> print(points.shape) (2, 3) >>> print(t0.indecies(points)) [14 3] >>> print(t1.indecies(points)) [58 15] >>> print(t2.indecies(points)) [235 62] ''' if points.ndim == 1: return self.index1(points) # Ok, delegate!! logger = self.logger # logger.setLevel(logging.INFO) if self.parent is None: return self._search_indecies(points) parent_index = self.parent.indecies(points) # (N,) array. Index of parent triangle for each point. logger.debug('index of parent = {}'.format(parent_index)) parent_index_set = set(parent_index) # A unique set for index of parent triangle index = np.zeros(points.shape[0], dtype=int) + len(self.tess_tris) # Here, points is assuemd 2D! Initial index is more than the allowed. for pindex in parent_index_set: # Instead of searching full parent triangles, only unique set is searched. index_of_points_for_pindex = np.where(parent_index == pindex)[0] # Index of points where parent index is pindex possible_indexes = self.index_map[pindex] # For each parent triangle, the derived indexes of triangles will be examined. t0, t1, t2, tc = self.tess_tris[possible_indexes] t0arr = issamehemisphere(points[index_of_points_for_pindex], t0[0], t0[1], t0[2]) tarr = (np.zeros_like(index, dtype=bool)) tarr[index_of_points_for_pindex] = t0arr index = np.where(tarr, possible_indexes[0], index) t1arr = issamehemisphere(points[index_of_points_for_pindex], t1[0], t1[1], t1[2]) tarr = (np.zeros_like(index, dtype=bool)) tarr[index_of_points_for_pindex] = t1arr index = np.where(tarr, possible_indexes[1], index) t2arr = issamehemisphere(points[index_of_points_for_pindex], t2[0], t2[1], t2[2]) tarr = (np.zeros_like(index, dtype=bool)) tarr[index_of_points_for_pindex] = t2arr index = np.where(tarr, possible_indexes[2], index) t3arr = np.logical_not(np.logical_or(np.logical_or(t0arr, t1arr), t2arr)) tarr = (np.zeros_like(index, dtype=bool)) tarr[index_of_points_for_pindex] = t3arr index = np.where(tarr, possible_indexes[3], index) return index
def _createLinkIcosahedron(self): ''' Create information table of triangle link for ndim=0. ''' ic = icosahedron.RegularIcosahedron() self.linklist = ic.getTriangleLinks() # tris is a tuple of Triangle instances.
[docs] def index1(self, point): ''' Return the index of the triangle that contains points :param point: (3,)-shaped numpy array that gives a point coordinate. :returns: Index of triangle that contains the ``point``. To confirm the functionality, this method is implemented. So that this is rather a development version, though, it can be used also for practical use. It is limited to 1 point evaluation. Use :meth:`indecies` for practical use. >>> lo = np.deg2rad(18.7) >>> la = np.deg2rad(-24.9) >>> point = np.array([np.cos(lo) * np.cos(la), np.sin(lo) * np.cos(la), np.sin(la)]) >>> t0 = RecurrentTessellatedUnitSphere(ndim=0) >>> print(t0.index1(point)) 14 >>> t1 = RecurrentTessellatedUnitSphere(ndim=1) >>> print(t1.index1(point)) 58 >>> t2 = RecurrentTessellatedUnitSphere(ndim=2) >>> print(t2.index1(point)) 235 Another point for testing. >>> lo = np.deg2rad(-78.1) >>> la = np.deg2rad(58.2) >>> point = np.array([np.cos(lo) * np.cos(la), np.sin(lo) * np.cos(la), np.sin(la)]) >>> t0 = RecurrentTessellatedUnitSphere(ndim=0) >>> print(t0.index1(point)) 3 >>> t1 = RecurrentTessellatedUnitSphere(ndim=1) >>> print(t1.index1(point)) 15 >>> t2 = RecurrentTessellatedUnitSphere(ndim=2) >>> print(t2.index1(point)) 62 ''' logger = self.logger # logger.setLevel(logging.INFO) logger.debug('index1 called with ndim={}'.format(self.ndim)) if self.parent is None: return self._search_index1(point) parent_index = self.parent.index1(point) logger.debug('index of parent = {}'.format(parent_index)) possible_indexes = self.index_map[parent_index] logger.debug('indecies of this = {}'.format(possible_indexes)) t0, t1, t2, tc = [self.tess_tris[i] for i in possible_indexes] # Four derived triangles assumed. logger.debug(' t0: {}'.format(t0)) logger.debug(' t1: {}'.format(t1)) logger.debug(' t2: {}'.format(t2)) logger.debug(' tc: {}'.format(tc)) # Here, obtain the index of triangle where the given point (piont) is contained. if issamehemisphere(point, t0[0], t0[1], t0[2]): ### t0 is the triangle that contains the first vertex of parent triangle. ### t0[0] is the vertex that forms the original triangle. ### t0[1] and t0[2] are the vertecies that are generated newly. ### So here it checks if the given point (point) is in the same side ### as t0[0] relative to the great circle formed by t0[1] and t0[2]. ### If yes, the given point is in the first triangle. ### If no, it is not. logger.debug('triangle: {}'.format(t0)) return possible_indexes[0] elif issamehemisphere(point, t1[0], t1[1], t1[2]): logger.debug('triangle: {}'.format(t1)) return possible_indexes[1] elif issamehemisphere(point, t2[0], t2[1], t2[2]): logger.debug('triangle: {}'.format(t2)) return possible_indexes[2] else: logger.debug('triangle: {}'.format(tc)) return possible_indexes[3]
def _new_triangles(self): ''' Create a tesselated sphere ''' ic = icosahedron.RegularIcosahedron() tris = ic.getNpTriangles() # tris is a tuple of Triangle instances. return tris def _split_triangles(self, tess_tris): new_tris = [] index_list = [] for tri in tess_tris: new_tess_tris = mesh_triangle(tri, output_normalize_factor=1.0) ### input_normalize_factor is not given for performance, ### and here assumed I the input vertexes have length=1. indices = tuple(range(len(new_tris), len(new_tris) + len(new_tess_tris))) new_tris += new_tess_tris index_list.append(indices) return new_tris, index_list
[docs] def indices_lonlat(self, lon_arr, lat_arr): ''' A wrapper of calling :meth:`indices`. :param lon_arr: ``np.array`` of longitude in degrees. Any angle acceptable. :param lat_arr: ``np.array`` of latitude in degrees. It should be -90 and 90 in theory, but no range checks are done. :returns: ``np.array`` of indecies each point belongs. >>> lo = np.array(18.7) >>> la = np.array(-24.9) >>> t0 = RecurrentTessellatedUnitSphere(ndim=0) >>> print(t0.indices_lonlat(lo, la)) 14 >>> t1 = RecurrentTessellatedUnitSphere(ndim=1) >>> print(t1.indices_lonlat(lo, la)) 58 >>> t2 = RecurrentTessellatedUnitSphere(ndim=2) >>> print(t2.indices_lonlat(lo, la)) 235 >>> lo = np.array([18.7, -78.1]) >>> la = np.array([-24.9, 58.2]) >>> print(t0.indices_lonlat(lo, la)) [14 3] >>> print(t1.indices_lonlat(lo, la)) [58 15] >>> print(t2.indices_lonlat(lo, la)) [235 62] ''' lon = np.deg2rad(lon_arr) lat = np.deg2rad(lat_arr) x = np.cos(lon) * np.cos(lat) y = np.sin(lon) * np.cos(lat) z = np.sin(lat) pnt = np.array([x, y, z]).T return self.indecies(pnt)
[docs] def search_area(self, bool_function, *args, strategy="array_full", **kwds): ''' Search the area that satisfies the given condition This is an interface providing area searches. Here area search means a search of an area/areas that satisfies the ``bool_function``. Several implementations are possible, depending on the bool_function types. Users can specify the implementation by ``strategy`` keyword. :param bool_function: Function that eats arrays and/or scalars (depending on strategy) of longitude and latitude returning an array or scalar of bool. Here longitude and latitude should be in degrees, and longitude range are -180 and 180 degrees. :keyword strategy: Implementation of strategy. "array_full", "array_link", "single_link", and "single_full" is supported as of now. Other keywords are forwarded to each implementation. Strategies: * "array_full": ``bool_function`` should eat ``np.array`` of longitude and latitude, returning the same shaped ``np.array``. All the triangles are examined. * "single_link": ``bool_function`` eats two argument, the longitude and latitude. By specifying "pre-known" coordinate that gives ``True`` in ``bool_function``, evaluations of near-by triangles are executed until all the possible neighbors become False. Thus, only one region can be identified. * "array_link": Similar to "single_link", though the ``bool_function`` can eat arrays of longitude and latitude parameters. * "single_full" The worst performance version. All the triangle is iterated by calling ``bool_function`` each by each. How to select strategies: To be written. ''' if strategy == 'array_full': return self.search_area_array_full(bool_function) elif strategy == 'array_link': return self.search_area_array_link(bool_function, args[0], **kwds) elif strategy == 'single_full': return self.search_area_single_full(bool_function) elif strategy == 'single_link': return self.search_area_single_link(bool_function, args[0], **kwds) else: raise ValueError('Strategy {} is unknown'.format(strategy))
[docs] def search_area_array_full(self, bool_function): ''' Search the area that satisfies the given condition This method returns a tuple of indecies of which the center of gravity satisfies the bool_function. The implementation can be called from the interface of :meth:`search_area` with a strategy of "array_full". In this implementation, The user have to provide a function that eats the arrays of longitude and latitude returning the boolean array. The calculation is done at once using numpy array, so that this implementaion is in general the fastest and most comprehensive. Drawback is that the implementation of bool function is a bit more complicated that the scaler arguments version (:meth:`search_area_array_full`). For example, if you want to get the area where north latitude is more than 66.7 degrees and positive longitude (0<=longitude<=180), you may imprement a bool function like >>> def in_east_arctic(lon, lat): ... return np.logical_and(lat > 66.7, lon >= 0) >>> sphere = RecurrentTessellatedUnitSphere(ndim=2) >>> idx_array = sphere.search_area_array_full(in_east_arctic) >>> print(idx_array) (0, 3, 16, 19, 32, 35) or, using general interface >>> idx_array = sphere.search_area(in_east_arctic, strategy='array_full') >>> print(idx_array) (0, 3, 16, 19, 32, 35) ''' cogs = self.cogs lons = np.rad2deg(np.arctan2(cogs[:, 1], cogs[:, 0])) lats = np.rad2deg(np.arcsin(cogs[:, 2])) bool_array = bool_function(lons, lats) idx = np.where(bool_array)[0] return tuple(sorted(idx))
[docs] def search_area_single_full(self, bool_function): ''' Search the area that satisfies the given condition This method returns a tuple of indecies of which the center of gravity satisfies the bool_function. The implementation can be called from the interface of :meth:`search_area` with a strategy of "single_full". In this implementation, the user have to provide a function that eats longitude (-180 to 180) and latitude scalars returning the boolean array. The evaluation is done over all traiangles each by each, so that the execution time will be very long. However, this is the most comprehensive solution, if the implementation of array version of bool function is hard to implement. For example, if you want to get the area where north latitude is more than 66.7 degrees and positive longitude (0<=longitude<=180), you may imprement a bool function like >>> def in_east_arctic(lon, lat): ... return (lat > 66.7) & (lon >= 0) # Scalar version >>> sphere = RecurrentTessellatedUnitSphere(ndim=2) >>> idx_array = sphere.search_area_single_full(in_east_arctic) >>> print(idx_array) (0, 3, 16, 19, 32, 35) or, using general interface >>> idx_array = sphere.search_area(in_east_arctic, strategy='single_full') >>> print(idx_array) (0, 3, 16, 19, 32, 35) ''' cogs = self.cogs lons = np.rad2deg(np.arctan2(cogs[:, 1], cogs[:, 0])) lats = np.rad2deg(np.arcsin(cogs[:, 2])) idx = [] for i in range(len(lons)): if bool_function(lons[i], lats[i]): idx.append(i) return tuple(sorted(idx))
[docs]def mesh_triangle(npTriangle, input_normalize_factor=None, output_normalize_factor=None): r''' A given triangle is meshed into four. A given triangle is separated into four triangles. :param npTriangle: NpTriangle object A triangle to be meshed. :keyword input_normalize_factor: If not ``None``, the input triangle's vertexes are normalized to the number given. :keyword output_normalize_factor: If not ``None``, the generated triangle's vertexes are normalized to the number given. :type output_normalize_factor: ``float`` or ``None``. :Returns: t0, t1, t2, t3 : New triangles Separated triangles. The separation is performed exactly following this diagram:: v0 ---- v5 ---- v2 \ / \ / v3 --- v4 \ / v1 Triangle (v0, v1, v2) is separated to (v0, v3, v5) (v1, v4, v3) (v2, v5, v4) (v3, v4, v5) All triangles will be co-planerly, unless ``output_normalize_factor`` is specified. >>> p0 = np.array([1., 0., 0.]) >>> p1 = np.array([0., 1., 0.]) >>> p2 = np.array([0., 0., 1.]) >>> torg = triangle.NpTriangle(p0, p1, p2) >>> t0, t1, t2, t3 = mesh_triangle(torg) >>> print(t0) # doctest: +NORMALIZE_WHITESPACE NpTriangle([1. 0. 0.] [0.5 0.5 0. ] [0.5 0. 0.5]) >>> print(t3) # doctest: +NORMALIZE_WHITESPACE NpTriangle([0.5 0.5 0. ] [0. 0.5 0.5] [0.5 0. 0.5]) If ``output_normalize_factor`` is given, all the points are normalized. >>> t0, t1, t2, t3 = mesh_triangle(torg, output_normalize_factor=1.) >>> print(t0) # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE NpTriangle([1. 0. 0.] [0.707... 0.707... 0. ] [0.707... 0. 0.707...]) ''' v0, v1, v2 = npTriangle.getVertices() if input_normalize_factor is not None: v0 = v0 / np.linalg.norm(v0) * input_normalize_factor v1 = v1 / np.linalg.norm(v1) * input_normalize_factor v2 = v2 / np.linalg.norm(v2) * input_normalize_factor v3 = (v0 + v1) / 2. v4 = (v1 + v2) / 2. v5 = (v2 + v0) / 2. if output_normalize_factor is not None: v3 /= np.linalg.norm(v3) * output_normalize_factor v4 /= np.linalg.norm(v4) * output_normalize_factor v5 /= np.linalg.norm(v5) * output_normalize_factor return (triangle.NpTriangle(v0, v3, v5), triangle.NpTriangle(v1, v4, v3), triangle.NpTriangle(v2, v5, v4), triangle.NpTriangle(v3, v4, v5))
[docs]def issamehemisphere(points, reference, split0, split1): r''' Check if the given points and refernce is in the same hemishere. :param points: Points to examine. (3,) np.array or (N, 3) np.array :param reference: Reference point :param split0: A point on the great circle of interest :param split1: Another point on the great circle of interest :returns: Scalar or (N,) shaped array of boolean. Consider a great circle defined by points ``split0`` and ``split1``. The great circle separates a sphere into two. If the given points to examine and the reference point is same hemisphere, True is returned. It is a user's respoinsibility to normalize all the points (for example to the unit vector). For example, you can consider the equator. >>> split0 = np.array((1., 0, 0)) >>> split1 = np.array((0, 1., 0)) >>> reference = np.array((0, 0, 1.)) For the above setting, the northern hemisphere (z>0) is True, and the southern hemisphere (z<0) is False. >>> print(issamehemisphere(np.array([0.5, 0.2, 0.3]), reference, split0, split1)) True >>> print(issamehemisphere(np.array([0.5, 0.2, -0.3]), reference, split0, split1)) False Multiple points can be given. >>> p0 = [0.3, 0.9, 0.2] >>> p1 = [0.3, -0.7, -0.5] >>> p = np.array([p0, p1]) >>> print(p.shape) (2, 3) >>> print(issamehemisphere(p, reference, split0, split1)) [ True False] Algorithm is as follows. We write the split points to be :math:`S_0, S_1` and the reference point to be :math:`R_0`. The vector of the point to examine to be :math:`P`. A normal vector of the great circle, :math:`N` is given by :math:`N=S_0\times S_1`. The sign of the inner product to :math:`N` will give the side of the vector of interest. Thus, the sign of :math:`(N\cdot P)\times(N\cdot R_0)` will give if :math:`P` and :math:`R_0` is in the same side (>0) or the other side (<0). ''' normvec = np.cross(split0, split1) if (normvec == [0, 0, 0]).all(): raise ValueError('A great circle cannot be defined.') curval = np.dot(points, normvec) refval = np.dot(reference, normvec) return curval * refval >= 0
import unittest import doctest
[docs]def doctests(): return unittest.TestSuite(( doctest.DocTestSuite(), ))
if __name__ == '__main__': unittest.main(defaultTest='doctests')