/
utils.py
142 lines (123 loc) · 4.11 KB
/
utils.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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
"""
Utils for making maps with ipyvolume.
This contains functions to convert between spherical latitude/longitude
and Cartesian x, y, z coordinates, and some functions to extract data
from GeoJSON (or TopoJSON, a slightly more compact version of GeoJSON).
"""
from typing import Tuple
from math import radians, log, sin, cos, tan, asin, atan2, pi
pi_div_180 = pi / 180
def latlon2xyz(
lat: float,
lon: float,
radius: float=1,
unit: str='deg'
) -> Tuple[float, float, float]:
"Convert lat/lon pair to Cartesian x/y/z triple."
if unit == 'deg':
lat = lat * pi_div_180
lon = lon * pi_div_180
cos_lat = cos(lat)
x = radius * cos_lat * sin(lon)
y = radius * sin(lat)
z = radius * cos_lat * cos(lon)
return (x, y, z)
def lonlat2xyz(lon, lat, radius=1, unit='deg'):
"Convert lat/lon pair to Cartesian x/y/z triple."
if unit == 'deg':
lat = lat * pi_div_180
lon = lon * pi_div_180
cos_lat = cos(lat)
x = radius * cos_lat * sin(lon)
y = radius * sin(lat)
z = radius * cos_lat * cos(lon)
return (x, y, z)
def xyz2latlon(
x: float,
y: float,
z: float,
radius: float=1,
unit: str='deg'
) -> Tuple[float, float]:
"Convert Cartesian x/y/z triple to lat/lon pair."
lat = asin(z / radius)
lon = atan2(y, x)
if unit == 'deg':
lat = lat / pi_div_180
lon = lon / pi_div_180
return (lat, lon)
def extract_coords(gj: dict):
"Yield all points in some GeoJSON/TopoJSON as [lon, lat] coordinate pairs."
gj_type = gj['type']
if gj_type == 'Point':
yield gj['coordinates']
elif gj_type in ['MultiPoint', 'LineString']:
for coord in gj_coords: # FIXME
yield coord
elif gj_type in ['MultiLineString', 'Polygon']:
for line in gj['coordinates']:
for coord in line:
yield coord
elif gj_type == 'MultiPolygon':
for poly in gj['coordinates']:
for line in poly:
for coord in line:
yield coord
elif gj_type == 'GeometryCollection':
for geom in gj['geometries']:
for coord in extract_coords(geom):
yield coord
elif gj_type == 'FeatureCollection':
for feat in gj['features']:
for coord in extract_coords(feat):
yield coord
elif gj_type == 'Feature':
geom = gj['geometry']
for coord in extract_coords(geom):
yield coord
else:
msg = f'Unkown GeoJSON/TopoJSON type: {gj_type}'
raise ValueError(msg)
def extract_lines(gj: dict):
"""
Yield lines in some GeoJSON/TopoJSON as lists of [lon, lat] coordinate pairs.
This ignores Point and MultiPoint because they don't form lines.
"""
gj_type = gj['type']
if gj_type in ['LineString']:
yield gj['coordinates']
elif gj_type in ['MultiLineString', 'Polygon']:
yield gj['coordinates']
elif gj_type == 'MultiPolygon':
for poly in gj['coordinates']:
for line in poly:
yield line
elif gj_type == 'GeometryCollection':
for geom in gj['geometries']:
for line in extract_lines(geom):
yield line
elif gj_type == 'FeatureCollection':
for feat in gj['features']:
for line in extract_lines(feat):
yield line
elif gj_type == 'Feature':
geom = gj['geometry']
for line in extract_lines(geom):
yield line
elif gj_type == 'Topology':
transform = gj.get('transform', {})
scale = transform.get('scale', [1.0, 1.0])
translate = transform.get('translate', [0.0, 0.0])
for arc in gj['arcs']:
line = []
prev = [0, 0]
for point in arc:
prev[0] += point[0]
prev[1] += point[1]
line.append((prev[0] * scale[0] + translate[0], prev[1] * scale[1] + translate[1]))
yield line
elif gj_type in ['Point', 'MultiPoint']:
pass
else:
msg = f'Unknown GeoJSON/TopoJSON type: {gj_type}'
raise ValueError(msg)