-
Notifications
You must be signed in to change notification settings - Fork 0
/
topocad.py
189 lines (162 loc) · 5.68 KB
/
topocad.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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
# -*- coding: utf-8 -*-
"""
topocad
"""
import cadquery as cq
import numpy as np
from math import radians, cos, sin, asin, sqrt
def haversine(coord1, coord2):
"""Calculate the great circle distance in kilometers between two points on
the earth (specified in decimal degrees)
Parameters
----------
coord1 : tuple
(lat, lon) for point1 in decimal degrees
coord2 : tuple
(lat, lon) for point2 in decimal degrees
Returns
-------
dist : float
Distance between points in km
"""
lat1, lon1 = coord1
lat2, lon2 = coord2
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles.
return c * r
def coarsen(arr, factor):
"""Coarsen a 2D topo array by a factor
Parameters
----------
arr : np.ndarray
2D (y, x) array of elevation values in meters (z)
factor : int
Factor by which to coarsen. If the arr is (100, 100) and factor=2, the
output will be (50, 50)
Returns
-------
out : np.ndarray
Coarsened array. Each pixel is the average of factor**2 pixels from the
input.
"""
idy = factor * (arr.shape[0] // factor)
idx = factor * (arr.shape[1] // factor)
out = arr[:idy, :idx].reshape((idy // factor, factor,
idx // factor, factor))
out = out.mean(axis=(1, 3))
return out
def get_xz(arr, idy, dx, subsample=1, x_scale=1, z_exag=1, z_adder=1):
"""Get a tuple of many (x,z) coordinates corresponding to one row (idy) in
the topo array.
Parameters
----------
arr : np.ndarray
2D (y, x) array of elevation values in meters (z)
idy : int
Axis=0 index (y-axis)
dx : float
Distance in km along the x-axis
subsample : int
Interval at which to sample the x-axis (reduce detail for lower compute
costs)
x_scale : float
How to scale the x-axis values in mm. e.g., if this 100, the
x-coordinates will range from (0, 100) mm
z_exag : float
Option to exaggerate the topography (multiplier). 1 makes everthing to
true scale
z_adder : float
This will be the minimum z-value and results in the minimum thickness
in mm of the work piece.
Returns
-------
xz : tuple
Tuple of (x, z) points for cadquery to spline. x values are longitude
points in x_scale units z values are elevation values that are
appropriately scaled to the x values and the input z_exag.
"""
x = np.linspace(0, x_scale, arr.shape[1])
dz = arr.max() - arr.min()
# dz / (dx*1000) = true_z_exag / x_scale # maintain aspect ratio
true_z_scale = x_scale * dz / (dx*1000)
z = (arr[idy] - arr.min()) / dz
z = z * true_z_scale * z_exag + z_adder
xz = tuple(zip(x, z))
xz = xz[::subsample]
return xz
def make_topo_cad(arr, dx, dy, x_scale, z_exag, z_adder, subsample=50,
spline=False):
"""Make a 3D CAD topography model
Parameters
----------
arr : np.ndarray
2D (y, x) array of elevation values in meters (z)
dx : float
Distance in km along the x-axis
dy : float
Distance in km along the y-axis
x_scale : float
How to scale the x-axis values in mm. e.g., if this 100, the
x-coordinates will range from (0, 100) mm
z_exag : float
Option to exaggerate the topography (multiplier). 1 makes everthing to
true scale
z_adder : float
This will be the minimum z-value and results in
the minimum thickness in mm of the work piece.
subsample : int
Interval at which to sample the x and y axes e.g.,
slice(None, None, subsample). A larger value reduces detail for lower
compute costs. Compute scales approximately linearly with the remaining
points in the array. Based on a simple test on an M1 macbook air, the
cost for spline=False is approx 8.34e-4 seconds/point,
cost for spline=True is approx a lot.
spline : bool
Connect points with a spline (True) or linearly (False). Spline is more
computationally expensive.
Returns
-------
wp : cadquery.Workplane
Cadquery Workplane object with the 3D model. Units are in mm and
dimensions are specified by the x_scale, z_exag, and z_adder inputs.
ppm : float
The model's points per mm in the x-axis (measure of resolution/detail)
aspect : str
The x:y aspect ratio of the model
"""
wp = cq.Workplane("XZ")
irows = range(0, len(arr), subsample)
y_scale = x_scale * dy / dx
y_offset = y_scale / len(irows)
for i, irow in enumerate(irows):
if i > 0:
wp = wp.workplane(offset=y_offset)
xz = get_xz(arr, irow, dx,
subsample=subsample,
x_scale=x_scale,
z_exag=z_exag,
z_adder=z_adder)
if spline:
wp = (wp
.lineTo(*xz[0])
.spline(xz[1:], includeCurrent=True)
.lineTo(xz[-1][0], 0)
.close()
)
else:
for ixz in xz:
wp = wp.lineTo(*ixz)
wp = wp.lineTo(xz[-1][0], 0).close()
if spline:
wp = wp.loft(combine=True, ruled=False)
else:
wp = wp.loft(combine=False, ruled=True)
ppm = arr.shape[1] / subsample / x_scale
aspect = f'{x_scale:.0f}:{y_scale:.0f}'
return wp, ppm, aspect