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

• RecurrentTessellatedUnitSphere: An implementation of the tessellated sphere, a spphere represented by multiple triangles with almost the same area size.

• 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 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?)

Code author: Yoshifumi Futaana

class irfpy.util.tessellatedsphere.TessellatedSphere(ndim=0, radius=1.0)

Bases: object

A tessellated sphere.

This class represents a tessellated sphere.

Warning

This class is for reference implementation using legacy irfpy.util functionality. Use 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.

The constractor.

Parameters:

• ndim (int) – The degree of tessellation of the sphere. Default is 0, the icosahedron.

• radius (float) – Radius of the sphere. Default is 1.

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

getTriangles() gives back all the triangles. Triangle is implemented as 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 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)

getTriangles()

Returns the tuple of the triangles that tessellate the sphere.

class irfpy.util.tessellatedsphere.RecurrentTessellatedUnitSphere(ndim=0)

Bases: object

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 \(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.

The constractor.

Parameters:

ndim (int) – The degree of tessellation of the sphere. Default is 0, the icosahedron.

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

getTriangles() gives back all the triangles. Triangle is implemented as 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]

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 call createLinkList().

tess_tris

Tessellated triangles (ndarray of triangle.NpTriangle)

Numpy object array of triangle.NpTriangle. (N,) shape, where N is the number of triangles.

getTriangles() is the better way to get access to this variable.

index_map

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 indecies() or indecies_lonlat().

The shape of this variable 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.

getTriangles()

Return the (reference to) list of tessellated triangles.

getTriangleVertexes()

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

getTriangleLonLat()

Return the tessellating triangle longitude-latitude.

Similar to 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

indecies(points)

Return the indecies of the triangle IDs that contain the given 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 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]

createLinkList()

Create information table of triangle link (neighboring triangles).

This method create the linklist attribute that includes the link information.

Note

What is the link and link list? Here I mean the link is the triangles that are sharing a same edge. In other words, triangles sharing a border is “linked”.

The link list is the list of this information.

For example, linklist [127] provides indexes of triangles that shares an edge to the triangle 127.

Therefore, the shape of the list is (N, 3), where N is the number of triangles.

>>> ts = RecurrentTessellatedUnitSphere(ndim=3)
>>> print(ts.getTriangleVertexes().shape)
(1280, 3, 3)
>>> ts.createLinkList()
>>> print(ts.linklist.shape)
(1280, 3)

The triangle #127 shares the edges with triangles #124, 125, and 126.

>>> print(ts.linklist[127])
[124 125 126]

The opposite shall be true, the triangle #124 for example should share an edge with #127

>>> print(ts.linklist[124])
[113 118 127]

index1(point)

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)

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)

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

search_area_array_full(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 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)

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 include 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 search_area_single_full() or search_area_single_link(). Though the number of evaluation is more than 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)

search_area_single_full(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 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)

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

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 or None.) – If not None, 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)

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

irfpy.util.tessellatedsphere.doctests()