-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathlongitude_grid_source.py
36 lines (33 loc) · 1.51 KB
/
longitude_grid_source.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# ïàðàëëåëè
import numpy as np
from vtk.util.numpy_support import numpy_to_vtk
from vtk.numpy_interface import algorithms as algs
from vtk.numpy_interface import dataset_adapter as dsa
# øèðîòû, íà êîòîðûõ îòîáðàæàþòñÿ ïàðàëëåëè (ìèí, ìàêñ, øàã)
lon_step_deg = 10
medirians = np.arange(-180, 180, lon_step_deg)
N_OF_MERIDIANS = medirians.shape[0]
RESOLUTION = 100 # êîëè÷åñòâî òî÷åê â êàæäîé îêðóæíîñòè ìåðèäèàíà
EARTH_RADIUS_KM = 6372 # ðàäèóñ ñôåðè÷åñêîé Çåìëè, êì + 1 êì
# ïîñòðîåíèå ïàðàëëåëåé
vtkpoints = vtk.vtkPoints()
meridians_steps = np.linspace(-80,80, RESOLUTION) # ãðàäóñíûå ìåðû òî÷åê â îäíîé ïàðàëëåëè
for idx_mer, lon_deg in enumerate(medirians):
# ïîñòðîåíèå îäíîé ïàðàëëåëè
lon_rad = radians(lon_deg) # øèðîòà ïàðàëëåëè â ðàäèàíàõ
# íàáîð ïàðàëëåëè èç òî÷åê
for idx_pt, lat_deg in enumerate(meridians_steps):
lat_rad = radians(lat_deg) # äîëãîòà â ðàäèàíàõ
x = EARTH_RADIUS_KM * cos(lat_rad) * cos(lon_rad)
y = EARTH_RADIUS_KM * cos(lat_rad) * sin(lon_rad)
z = EARTH_RADIUS_KM * sin(lat_rad)
vtkpoints.InsertPoint(idx_mer*RESOLUTION + idx_pt, x,y,z)
output.SetPoints(vtkpoints)
print("vtkpoints = ", vtkpoints)
output.Allocate(N_OF_MERIDIANS, 1)
for idx_mer, lon_deg in enumerate(medirians):
vtkpolyline = vtk.vtkPolyLine()
vtkpolyline.GetPointIds().SetNumberOfIds(RESOLUTION)
for i in range(0,RESOLUTION):
vtkpolyline.GetPointIds().SetId(i, i + RESOLUTION * idx_mer)
output.InsertNextCell(vtkpolyline.GetCellType(), vtkpolyline.GetPointIds())