irfpy.util.tessellatedsphere
¶
Tessellating a sphere by triangles.
Code author: Yoshifumi Futaana
- class irfpy.util.tessellatedsphere.TessellatedSphere(ndim=None, radius=1.0)[source]¶
Bases:
object
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.
- Parameters
ndim (
int
or None) – Dimension of the tessellated sphere. IfNone
is given, ndim=0 is used as default.radius (
float
) – Radius of the circumscribe sphere.
>>> 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
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)
Constructor.
- class irfpy.util.tessellatedsphere.RecurrentTessellatedUnitSphere(ndim=0)[source]¶
Bases:
object
Tessellated sphere with keeping information on the parent-child relation
- logger = <Logger irfpy.util.tessellatedsphere.RecurrentTessellatedUnitSphere (DEBUG)>¶
- linklist¶
Contains information on triangle links.
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 callcreateLinkList()
.
- tess_tris¶
Tessellated triangles (ndarray of
triangle.NpTriangle
)Numpy object array of
triangle.NpTriangle
. (N,) shape, where N is the number of triangles.
- index_map¶
Index map from parent triangle to the self triangles.
The shape is (\(N_p\), the number of triangles of parent tessellation). If
ndim
is 0,None
is set.For example,
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.
- indecies(points)[source]¶
Return the indecies of the triangles that contain the points.
- Parameters
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
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]
- createLinkList()[source]¶
Create information table of triangle link.
This method create the
tri_link
attribute that includes the link information.For example,
tri_link
[127] is a tuple of indeces of triangles that shares an edge to the triangle 127.
- index1(point)[source]¶
Return the index of the triangle that contains points
- Parameters
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
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
- indices_lonlat(lon_arr, lat_arr)[source]¶
A wrapper of calling
indices()
.- Parameters
lon_arr –
np.array
of longitude in degrees. Any angle acceptable.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]
- search_area(bool_function, *args, strategy='array_full', **kwds)[source]¶
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.- Parameters
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.
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 eatnp.array
of longitude and latitude, returning the same shaped
np.array
. All the triangles are examined.
- “array_full”:
- “single_link”:
bool_function
eats two argument, the longitude and latitude. By specifying “pre-known” coordinate that gives
True
inbool_function
, evaluations of near-by triangles are executed until all the possible neighbors become False. Thus, only one region can be identified.
- “single_link”:
- “array_link”: Similar to “single_link”, though the
bool_function
can eat arrays of longitude and latitude parameters.
- “array_link”: Similar to “single_link”, though the
- “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.
- search_area_array_full(bool_function)[source]¶
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
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 (
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)
- search_area_array_link(bool_function, start_lonlat, maxeval=- 1)[source]¶
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
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
start_lonlat
argument. If the triangle that includestart_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 bymaxeval
, a ValueError is raised. Ifmaxeval
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. Ifmaxeval
is reached, ValueError is raised.The benefit of this implementation is rather small amount of evaluation needed for evaluation compared to
search_area_single_full()
orsearch_area_single_link()
. Though the number of evaluation is more thansearch_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)
- search_area_single_full(bool_function)[source]¶
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
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)
- search_area_single_link(bool_function, start_lonlat, maxeval=- 1)[source]¶
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
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
start_lonlat
argument. If the triangle that includestart_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 bymaxeval
, a ValueError is raised. Ifmaxeval
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. Ifmaxeval
is reached, ValueError is raised.The benefit of this implementation is rather small amount of evaluation needed for evaluation compared to
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)
- irfpy.util.tessellatedsphere.mesh_triangle(npTriangle, input_normalize_factor=None, output_normalize_factor=None)[source]¶
A given triangle is meshed into four.
A given triangle is separated into four triangles.
- Parameters
npTriangle – NpTriangle object A triangle to be meshed.
input_normalize_factor – If not
None
, the input triangle’s vertexes are normalized to the number given.output_normalize_factor (
float
orNone
.) – If notNone
, the generated triangle’s vertexes are normalized to the number given.
- Returns
- t0, t1, t2, t3New 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) NpTriangle([1. 0. 0.] [0.5 0.5 0. ] [0.5 0. 0.5]) >>> print(t3) 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) NpTriangle([1. 0. 0.] [0.707... 0.707... 0. ] [0.707... 0. 0.707...])
- irfpy.util.tessellatedsphere.issamehemisphere(points, reference, split0, split1)[source]¶
Check if the given points and refernce is in the same hemishere.
- Parameters
points – Points to examine. (3,) np.array or (N, 3) np.array
reference – Reference point
split0 – A point on the great circle of interest
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
andsplit1
. 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 \(S_0, S_1\) and the reference point to be \(R_0\). The vector of the point to examine to be \(P\). A normal vector of the great circle, \(N\) is given by \(N=S_0\times S_1\). The sign of the inner product to \(N\) will give the side of the vector of interest. Thus, the sign of \((N\cdot P)\times(N\cdot R_0)\) will give if \(P\) and \(R_0\) is in the same side (>0) or the other side (<0).