''' Tessellating a sphere by triangles.
.. 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:
''' Generate a tessellated sphere.
*ndim* is the number of division. If *ndim*\ =0, it is just a
regular icosahedron. If *ndim*\ =1, each triangle of the icosahedron
is separated into 4 by choosing the central point of the triangle
normalized to the sphere.
This means that the number of the triangles is 80.
Each increament of *ndim* makes one more division,
making 4x original triangles.
:keyword ndim: Dimension of the tessellated sphere.
If ``None`` is given, *ndim*\ =0 is used as default.
:type ndim: ``int`` or None
:keyword radius: Radius of the circumscribe sphere.
:type radius: ``float``
>>> ts = TessellatedSphere(ndim=0)
>>> tris = ts.getTriangles()
>>> print(len(tris))
20
>>> 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
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)
'''
def __init__(self, ndim=None, radius=1.0):
''' Constructor.
'''
if ndim is None:
self.ndim = 0
else:
self.ndim = ndim
self.radius = radius
self.tris = self.makeTessellation(self.ndim)
[docs] def makeTessellation(self, ndim):
''' Make *ndim*\ -order tessellation.
:param ndim: Dimension to be tessellated.
: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):
''' Given list of the triangles are separated.
All the trinagles in the given list is separated into 4, and
the separated triangles are returned as a list.
Normally one do not need to call this method.
'''
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):
''' A given triangle is separated 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
Separated triangles.
The separation is like this ::
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)
Here, v3, v4, and v5 is normalized.
>>> 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
Normally one do not need to call this method.
'''
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 reference to the tuple of triangles.
'''
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:
''' Tessellated sphere with keeping information on the parent-child relation
'''
logger = logging.getLogger(__name__ + '.RecurrentTessellatedUnitSphere')
def __init__(self, ndim=0):
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.
'''
self.index_map = None
''' Index map from parent triangle to the self triangles.
The shape 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))
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 triangles that contain the 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 ``points``.
For single point sample (OBS: this is delegated to :meth:`index1`),
>>> 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=np.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=np.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=np.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=np.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=np.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 createLinkList(self):
''' Create information table of triangle link.
This method create the :attr:`tri_link` attribute
that includes the link information.
For example, :attr:`tri_link` [127] is a tuple of
indeces of triangles that shares an edge to the
triangle 127.
'''
if self.ndim == 0:
self._createLinkIcosahedron()
return
else:
self.parent.createLinkList()
parent_linklist = self.parent.linklist
n = len(self.tess_tris)
self.linklist = np.zeros([n, 3]) + n # Initialize the list with "unrealistic" value.
n_p = len(self.parent.tess_tris)
# Probably very slow way, but aim to work.
for i_parent in range(n_p):
derived_triangle_indecies = self.index_map[i_parent]
### For the last one is simple, because this is the center!
self.linklist[derived_triangle_indecies[3], :] = derived_triangle_indecies[0:3]
### Parent link is important!
parent_link = parent_linklist[i_parent] # You get (3,) array of indecies
### Corresponding child link
possible_derived_link = self.index_map[parent_link][:, :3].flatten() # (9,) indecies, possible triangles that are neighboring.
possible_neighbor_triangles = self.tess_tris[possible_derived_link]
triangle_coords = [tri.getVertices() for tri in possible_neighbor_triangles] # (9, 3, 3) array
# print(triangle_coords)
for i, ti in enumerate(derived_triangle_indecies[:3]):
mytriangle = self.tess_tris[ti]
v0, v1, v2 = mytriangle.getVertices()
# print(mytriangle.getVertices())
sharelist = []
for vj in (v0, v1, v2):
triangle_coords_vj = triangle_coords - vj # (9, 3, 3) array
triangle_coords_norm_vj = np.linalg.norm(triangle_coords_vj, axis=2) # (9, 3) array
share_vj = (triangle_coords_norm_vj == 0.0) # (9, 3) booelan array
share_vj = share_vj.any(axis=1) # (9,) array
# if len(np.where(share_v0)[0]) != 2:
# raise RuntimeError('Check?')
sharelist.append(share_vj)
sharelist = np.array(sharelist) # (3, 9) boolean array
# print(sharelist)
sharelist = np.where(sharelist, 1, 0)
sharelist = sharelist.sum(axis=0)
linktriindex = np.where(sharelist == 2)[0]
# print(linktriindex)
self.linklist[derived_triangle_indecies[i], :2] = possible_derived_link[linktriindex] # Neighboring in other parent triangles
self.linklist[derived_triangle_indecies[i], 2] = derived_triangle_indecies[3] # Center in the same parent riangle
self.linklist = np.array(self.linklist, dtype=np.int64)
[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_array_link(self, bool_function, start_lonlat, maxeval=-1):
''' Search the area that satisfies the given condition
.. note::
Beta version. Not very well tested. Use carefully.
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_link".
In this implementation,
the user have to provide a function that eats
numpy arrays of longitude (-180 to 180) and latitude
returning the boolean array.
The evaluation is done in two steps.
In the first step, the method try to find a triangle
that satisfies the given condition with a help of
:attr:`start_lonlat` argument.
If the triangle that include :attr:`start_lonlat`
satisfies the given condition, it is considered as
a start triangle.
If not, neiboring triangles are evaluated until the
given condition is met.
If the evaluation of triangles exceeds the number given by
``maxeval``, a ValueError is raised.
If ``maxeval`` is -1, the evaluations will be done until all
the triangles are evaluated.
When the start triangle is found, the second step starts.
The neighboring triangles are evaluated recursively
until all the neighboring triangles unsatisfy the condition.
``maxeval`` parameters may also be used in the second step.
If ``maxeval`` is reached, ValueError is raised.
The benefit of this implementation is rather small amount of
evaluation needed for evaluation compared to
:meth:`search_area_single_full` or :meth:`search_area_single_link`.
Though the number of evaluation is more than
:meth:`search_area_array_full`,
each evaluation is a small chunk of arrays.
Thus, the total exec time will be smaller, in particular
for larger ndims.
However, limitation is that the search can find only one
connected area, and still the execution time can vary strongly
depending on the conditions etc.
Here is an example, if you want to get the area of the arctic
and antarctic circles, namely, the latitude is >66.7 or <-66.7,
and positive longitude (0<=longitude<=180),
you may imprement a bool function like
>>> def in_east_arctic(lon, lat):
... return np.logical_and(np.abs(lat) > 66.7, lon >= 0)
>>> sphere = RecurrentTessellatedUnitSphere(ndim=2)
>>> sphere.createLinkList()
>>> idx_array = sphere.search_area_array_link(in_east_arctic, [130., 75.4])
>>> print(idx_array)
(0, 3, 16, 19, 32, 35)
or, using general interface
>>> idx_array = sphere.search_area(in_east_arctic, [130., 75.4], strategy='array_link', maxeval=100000)
>>> print(idx_array)
(0, 3, 16, 19, 32, 35)
Note that even you want to look at the southern hemisphere,
it cannot be found. This is because the method look for only
the nearest region of ``start_lonlat``.
>>> idx_array = sphere.search_area(in_east_arctic, [120., -75.4], strategy='array_link', maxeval=100000)
>>> print(idx_array)
(244, 247, 260, 263)
'''
if self.linklist is None:
raise ValueError('Link list not available. Use #createLinkList() first.')
logger = self.logger
cogs = self.cogs
lons = np.rad2deg(np.arctan2(cogs[:, 1], cogs[:, 0]))
lats = np.rad2deg(np.arcsin(cogs[:, 2]))
lon, lat = np.deg2rad(start_lonlat)
p = np.array([np.cos(lon) * np.cos(lat),
np.sin(lon) * np.cos(lat),
np.sin(lat)])
stoplen = len(self.tess_tris) - maxeval
if maxeval < 0:
stoplen = 0
### First step
idx0 = self.index1(p)
set_unknown = set(range(len(self.tess_tris)))
set_true = set([])
set_false = set([])
set_candidate = set([idx0, ])
set_unknown.remove(idx0)
# print_status = lambda: logger.debug('''UNKNOWN: {}
#TRUE: {}
#FALSE: {}
#CANDIDATE: {}'''.format(set_unknown, set_true, set_false, set_candidate))
print_status = lambda: None
logger.debug('Start Step 1 with index {}'.format(idx0))
while len(set_unknown) + len(set_candidate) > stoplen:
print_status()
idxs = list(set_candidate)
result_of_candidate = bool_function(lons[idxs], lats[idxs])
# logger.debug(' EXAM: {}'.format(result_of_candidate))
if result_of_candidate.any(): # If any solution is found.
logger.debug(' Meet!')
idx_true = set(np.array(idxs)[np.where(result_of_candidate)[0]])
idx_false = set_candidate - idx_true
set_true |= idx_true
set_false |= idx_false
possible_links = []
for i in idx_true:
link = self.linklist[i]
possible_links += link.tolist()
set_candidate = set(possible_links)
set_candidate -= set_true
set_candidate -= set_false
set_unknown -= set_candidate
print_status()
break
else: # No solution is found.
logger.debug(' Unmeet!')
possible_links = []
for i in set_candidate:
link = self.linklist[i]
possible_links += link.tolist()
set_false |= set_candidate
set_candidate = set(possible_links)
print_status()
else:
raise ValueError('Maxeval reached (1st step)')
### Second step
logger.debug('Start Step 2')
while len(set_unknown) + len(set_candidate) > stoplen:
print_status()
if len(set_candidate) == 0:
break
idxs = list(set_candidate)
result_of_candidate = bool_function(lons[idxs], lats[idxs])
# logger.debug(' EXAM: {}'.format(result_of_candidate))
if result_of_candidate.any(): # If any of them is true
logger.debug(' Meet!')
idx_true = set(np.array(idxs)[np.where(result_of_candidate)[0]])
idx_false = set_candidate - idx_true
set_true |= idx_true
set_false |= idx_false
possible_links = []
for i in set_true:
link = self.linklist[i]
possible_links += link.tolist()
set_candidate = set(possible_links)
set_candidate -= set_true
set_candidate -= set_false
set_unknown -= set_candidate
print_status()
else:
logger.debug(' Unmeet!')
set_candidate = {}
print_status()
break # No other candidate!
else:
raise ValueError('Maxeval reached (2nd step)')
print_status()
return tuple(sorted(set_true))
[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 search_area_single_link(self, bool_function, start_lonlat, maxeval=-1):
''' Search the area that satisfies the given condition
.. note::
Beta version. Not very well tested. Use carefully.
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_link".
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 in two steps.
In the first step, the method try to find a triangle
that satisfies the given condition with a help of
:attr:`start_lonlat` argument.
If the triangle that include :attr:`start_lonlat`
satisfies the given condition, it is considered as
a start triangle.
If not, neiboring triangles are evaluated until the
given condition is met.
If the evaluation of triangles exceeds the number given by
``maxeval``, a ValueError is raised.
If ``maxeval`` is -1, the evaluations will be done until all
the triangles are evaluated.
When the start triangle is found, the second step starts.
The neighboring triangles are evaluated recursively
until all the neighboring triangles unsatisfy the condition.
``maxeval`` parameters may also be used in the second step.
If ``maxeval`` is reached, ValueError is raised.
The benefit of this implementation is rather small amount of
evaluation needed for evaluation compared to
:meth:`search_area_single_full`.
However, limitation is that the search can find only one
connected area, and still the execution time can vary strongly
depending on the conditions etc.
Here is an example, if you want to get the area of the arctic
and antarctic circles, namely, the latitude is >66.7 or <66.7,
and positive longitude (0<=longitude<=180),
you may imprement a bool function like
>>> def in_east_arctic(lon, lat):
... return (abs(lat) > 66.7) & (lon >= 0) # Scalar version
>>> sphere = RecurrentTessellatedUnitSphere(ndim=2)
>>> sphere.createLinkList()
>>> idx_array = sphere.search_area_single_link(in_east_arctic, [130., 75.4])
>>> print(idx_array)
(0, 3, 16, 19, 32, 35)
or, using general interface
>>> idx_array = sphere.search_area(in_east_arctic, [130., 75.4], strategy='single_link', maxeval=100000)
>>> print(idx_array)
(0, 3, 16, 19, 32, 35)
Note that even you want to look at the southern hemisphere,
it cannot be found. This is because the method look for only
the nearest region of ``start_lonlat``.
>>> idx_array = sphere.search_area(in_east_arctic, [120., -75.4], strategy='single_link', maxeval=100000)
>>> print(idx_array)
(244, 247, 260, 263)
'''
if self.linklist is None:
raise ValueError('Link list not available. Use #createLinkList() first.')
logger = self.logger
# logger.setLevel(logging.WARN)
cogs = self.cogs
lons = np.rad2deg(np.arctan2(cogs[:, 1], cogs[:, 0]))
lats = np.rad2deg(np.arcsin(cogs[:, 2]))
lon, lat = np.deg2rad(start_lonlat)
p = np.array([np.cos(lon) * np.cos(lat),
np.sin(lon) * np.cos(lat),
np.sin(lat)])
stoplen = len(self.tess_tris) - maxeval
if maxeval < 0:
stoplen = 0
### First step
idx0 = self.index1(p)
set_unknown = set(range(len(self.tess_tris)))
set_true = set([])
set_false = set([])
set_candidate = set([idx0, ])
set_unknown.remove(idx0)
logger.debug('Start Step 1')
while len(set_unknown) + len(set_candidate) > stoplen:
idx = set_candidate.pop()
logger.debug('==Index {}=='.format(idx))
if bool_function(lons[idx], lats[idx]):
logger.debug(' Meet!')
set_unknown = set_unknown.union(set_candidate)
set_candidate = set([])
possible_links = self.linklist[idx]
for i in possible_links:
if i in set_unknown:
set_unknown.remove(i)
set_candidate.add(i)
set_true = set([idx, ])
break
else:
logger.debug(' Unmeet!')
set_false.add(idx)
possible_links = self.linklist[idx]
for i in possible_links:
if i in set_unknown:
set_unknown.remove(i)
set_candidate.add(i)
else:
raise ValueError('Maxeval reached (1st step)')
### Second step
logger.debug('Start Step 2')
while len(set_unknown) + len(set_candidate) > stoplen:
if len(set_candidate) == 0:
break
idx = set_candidate.pop()
logger.debug('==Index {}=='.format(idx))
if bool_function(lons[idx], lats[idx]):
logger.debug(' Meet!')
possible_links = self.linklist[idx]
for i in possible_links:
if i in set_unknown:
set_unknown.remove(i)
set_candidate.add(i)
set_true.add(idx)
else:
logger.debug(' Unmeet!')
set_false.add(idx)
else:
raise ValueError('Maxeval reached (2nd step)')
return tuple(sorted(set_true))
[docs]def mesh_triangle(npTriangle, input_normalize_factor=None, output_normalize_factor=None):
''' 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')