-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Attemp to add spherical indexing system
- Loading branch information
1 parent
f0852f8
commit e889f96
Showing
9 changed files
with
447 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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) |
Binary file not shown.
Binary file added
BIN
+3.25 KB
stemflow/utils/sphere/__pycache__/coordinate_transform.cpython-39.pyc
Binary file not shown.
Binary file not shown.
Binary file added
BIN
+3.01 KB
stemflow/utils/sphere/__pycache__/lonlat_spherical_transformer.cpython-39.pyc
Binary file not shown.
Oops, something went wrong.