Attemp to add spherical indexing system

.gitignore

@@ -15,3 +15,4 @@ stemflow/__pycache__/*
 *.docx
 */references/*
 .github/workflows/publish-pypi.yml
+**__pycache__**

stemflow/gridding/Q_blocks.py

@@ -1,5 +1,10 @@
+"""I call this Q_blocks because they are essential blocks for QTree methods"""
+
 from typing import Tuple, Union
 
+from ..utils.sphere.coordinate_transform import lonlat_spherical_transformer
+from ..utils.sphere.distance import spherical_distance_from_coordinates
+
 
 class Point:
     """A Point class for recording data points"""
@@ -14,7 +19,12 @@ class Node:
     """A tree-like division node class"""
 
     def __init__(
-        self, x0: Union[float, int], y0: Union[float, int], w: Union[float, int], h: Union[float, int], points: Point
+        self,
+        x0: Union[float, int],
+        y0: Union[float, int],
+        w: Union[float, int],
+        h: Union[float, int],
+        points: list[Point],
     ):
         self.x0 = x0
         self.y0 = y0
@@ -42,3 +52,46 @@ def __init__(self, x_index, y_index, x_range, y_range):
         self.x_range = x_range
         self.y_range = y_range
         self.points = []
+
+
+class Sphere_Point:
+    """A Point class for recording data points"""
+
+    def __init__(self, index, inclination, azimuth):
+        self.inclination = inclination
+        self.azimuth = azimuth
+        self.index = index
+
+
+class Sphere_Face:
+    """A tree-like division triangle node class for Sphere Quadtree"""
+
+    def __init__(
+        self,
+        x0: Union[float, int],
+        y0: Union[float, int],
+        azimuth1: Union[float, int],
+        inclination1: Union[float, int],
+        azimuth2: Union[float, int],
+        inclination2: Union[float, int],
+        azimuth3: Union[float, int],
+        inclination3: Union[float, int],
+        points: list[Sphere_Point],
+    ):
+        self.x0 = x0
+        self.y0 = y0
+
+        self.distance = spherical_distance_from_coordinates(
+            inclination1, azimuth1, inclination2, azimuth2, radius=6371.0
+        )
+        self.points = points
+        self.children = []
+
+    def get_width(self):
+        return self.width
+
+    def get_height(self):
+        return self.height
+
+    def get_points(self):
+        return self.points

stemflow/gridding/Sphere_QTree.py

@@ -0,0 +1,175 @@
+# import math
+# import os
+# import warnings
+# from collections.abc import Sequence
+
+# # from multiprocessing import Pool
+# from typing import Tuple, Union
+
+# import matplotlib.patches as patches
+# import matplotlib.pyplot as plt  # plotting libraries
+# import numpy as np
+# import pandas
+# import pandas as pd
+
+# from ..utils.generate_soft_colors import generate_soft_color
+# from ..utils.sphere.coordinate_transform import lonlat_spherical_transformer
+# from ..utils.sphere.distance import spherical_distance_from_coordinates
+# from ..utils.sphere.Icosahedron import get_earth_Icosahedron_vertices_and_faces
+# from ..utils.validation import check_random_state
+# from .Q_blocks import Grid, Node, Point
+
+# os.environ["MKL_NUM_THREADS"] = "1"
+# os.environ["NUMEXPR_NUM_THREADS"] = "1"
+# os.environ["OMP_NUM_THREADS"] = "1"
+
+# warnings.filterwarnings("ignore")
+
+
+# def recursive_subdivide(
+#     node: Node,
+#     grid_len_upper_threshold: Union[float, int],
+#     grid_len_lower_threshold: Union[float, int],
+#     points_lower_threshold: Union[float, int],
+# ):
+#     """recursively subdivide the grids
+
+#     Args:
+#         node:
+#             node input
+#         grid_len_lon_upper_threshold:
+#             force divide if grid longitude larger than the threshold
+#         grid_len_lon_lower_threshold:
+#             stop divide if grid longitude **will** be below than the threshold
+#         grid_len_lat_upper_threshold:
+#             force divide if grid latitude larger than the threshold
+#         grid_len_lat_lower_threshold:
+#             stop divide if grid latitude **will** be below than the threshold
+
+#     """
+#     raise NotImplementedError()
+
+#     if node.width / 2 < grid_len_lower_threshold:
+#         # The width and height will be the same. So only test one.
+#         return
+
+#     w_ = float(node.width / 2)
+
+#     p = contains(node.x0, node.y0, w_, h_, node.points)
+#     x1 = Node(node.x0, node.y0, w_, h_, p)
+#     recursive_subdivide(
+#         x1,
+#         grid_len_lon_upper_threshold,
+#         grid_len_lon_lower_threshold,
+#         grid_len_lat_upper_threshold,
+#         grid_len_lat_lower_threshold,
+#         points_lower_threshold,
+#     )
+
+#     p = contains(node.x0, node.y0 + h_, w_, h_, node.points)
+#     x2 = Node(node.x0, node.y0 + h_, w_, h_, p)
+#     recursive_subdivide(
+#         x2,
+#         grid_len_lon_upper_threshold,
+#         grid_len_lon_lower_threshold,
+#         grid_len_lat_upper_threshold,
+#         grid_len_lat_lower_threshold,
+#         points_lower_threshold,
+#     )
+
+#     p = contains(node.x0 + w_, node.y0, w_, h_, node.points)
+#     x3 = Node(node.x0 + w_, node.y0, w_, h_, p)
+#     recursive_subdivide(
+#         x3,
+#         grid_len_lon_upper_threshold,
+#         grid_len_lon_lower_threshold,
+#         grid_len_lat_upper_threshold,
+#         grid_len_lat_lower_threshold,
+#         points_lower_threshold,
+#     )
+
+#     p = contains(node.x0 + w_, node.y0 + h_, w_, h_, node.points)
+#     x4 = Node(node.x0 + w_, node.y0 + h_, w_, h_, p)
+#     recursive_subdivide(
+#         x4,
+#         grid_len_lon_upper_threshold,
+#         grid_len_lon_lower_threshold,
+#         grid_len_lat_upper_threshold,
+#         grid_len_lat_lower_threshold,
+#         points_lower_threshold,
+#     )
+
+#     for ch_node in [x1, x2, x3, x4]:
+#         if len(ch_node.points) <= points_lower_threshold:
+#             if not ((node.width > grid_len_lon_upper_threshold) or (node.height > grid_len_lat_upper_threshold)):
+#                 return
+
+#     node.children = [x1, x2, x3, x4]
+
+
+# def contains(x, y, w, h, points):
+#     """return list of points within the grid"""
+#     pts = []
+#     for point in points:
+#         if point.x >= x and point.x <= x + w and point.y >= y and point.y <= y + h:
+#             pts.append(point)
+#     return pts
+
+
+# def find_children(node):
+#     """return children nodes of this node"""
+#     if not node.children:
+#         return [node]
+#     else:
+#         children = []
+#         for child in node.children:
+#             children += find_children(child)
+#     return children
+
+
+# class Sphere_QTree:
+#     """A spherical Quadtree class"""
+
+#     def __init__(
+#         self,
+#         grid_len_upper_threshold: Union[float, int],
+#         grid_len_lower_threshold: Union[float, int],
+#         points_lower_threshold: int,
+#         lon_lat_equal_grid: bool = True,
+#         rotation_angle: Union[float, int] = 0,
+#         calibration_point_x_jitter: Union[float, int] = 0,
+#         calibration_point_y_jitter: Union[float, int] = 0,
+#     ):
+#         pass
+
+#     def add_lon_lat_data(self, indexes: Sequence, x_array: Sequence, y_array: Sequence):
+#         """Store input lng lat data and transform to **Point** object
+
+#         Parameters:
+#             indexes: Unique identifier for indexing the point.
+#             x_array: longitudinal values.
+#             y_array: latitudinal values.
+
+#         """
+#         if not len(x_array) == len(y_array) or not len(x_array) == len(indexes):
+#             raise ValueError("input longitude and latitude and indexes not in same length!")
+
+#         data = np.array([x_array, y_array]).T
+#         angle = self.rotation_angle
+#         r = angle / 360
+#         theta = r * np.pi * 2
+#         rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
+#         data = data @ rotation_matrix
+#         lon_new = (data[:, 0] + self.calibration_point_x_jitter).tolist()
+#         lat_new = (data[:, 1] + self.calibration_point_y_jitter).tolist()
+
+#         for index, lon, lat in zip(indexes, lon_new, lat_new):
+#             self.points.append(Point(index, lon, lat))
+
+#     def generate_gridding_params(self):
+#         """For completeness"""
+#         pass
+
+#     def get_points(self):
+#         """For completeness"""
+#         pass

stemflow/utils/sphere/Icosahedron.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+from .coordinate_transform import cartesian_3D_to_lonlat
+
+
+def get_Icosahedron_vertices():
+    phi = (1 + np.sqrt(5)) / 2
+    vertices = np.array(
+        [
+            (phi, 1, 0),
+            (phi, -1, 0),
+            (-phi, -1, 0),
+            (-phi, 1, 0),
+            (1, 0, phi),
+            (-1, 0, phi),
+            (-1, 0, -phi),
+            (1, 0, -phi),
+            (0, phi, 1),
+            (0, phi, -1),
+            (0, -phi, -1),
+            (0, -phi, 1),
+        ]
+    )
+    return vertices
+
+
+def calc_and_judge_distance(v1, v2, v3):
+    d1 = np.sum((np.array(v1) - np.array(v2)) ** 2) ** (1 / 2)
+    d2 = np.sum((np.array(v1) - np.array(v3)) ** 2) ** (1 / 2)
+    d3 = np.sum((np.array(v2) - np.array(v3)) ** 2) ** (1 / 2)
+    if d1 == d2 == d3 == 2:
+        return True
+    else:
+        return False
+
+
+def get_Icosahedron_faces():
+    vertices = get_Icosahedron_vertices()
+
+    face_list = []
+    for vertice1 in vertices:
+        for vertice2 in vertices:
+            for vertice3 in vertices:
+                same_face = calc_and_judge_distance(vertice1, vertice2, vertice3)
+                if same_face:
+                    the_face_set = set([tuple(vertice1), tuple(vertice2), tuple(vertice3)])
+                    if the_face_set not in face_list:
+                        face_list.append(the_face_set)
+
+    face_list = np.array([list(i) for i in face_list])
+    return face_list
+
+
+def get_earth_Icosahedron_vertices_and_faces():
+    # earth_radius_km=6371.0
+    # get Icosahedron vertices and faces
+    vertices = get_Icosahedron_vertices()
+    face_list = get_Icosahedron_faces()
+
+    # Scale: from 2 to 6371
+    # scale_ori = (np.sum(vertices**2, axis=1)**(1/2))[0]
+    # # Scale vertices and face_list to km
+    # vertices = vertices * (earth_radius_km/scale_ori)
+    # face_list = face_list * (earth_radius_km/scale_ori)
+
+    vertices_lng, vertices_lat = cartesian_3D_to_lonlat(vertices[:, 0], vertices[:, 1], vertices[:, 2])
+
+    faces_lng, faces_lat = cartesian_3D_to_lonlat(face_list[:, :, 0], face_list[:, :, 1], face_list[:, :, 2])
+
+    return np.stack([vertices_lng, vertices_lat], axis=-1), np.stack([faces_lng, faces_lat], axis=-1)