|
1
|
|
|
# Copyright (c) 2008-2016 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
""" |
|
5
|
|
|
============================= |
|
6
|
|
|
Natural Neighbor Verification |
|
7
|
|
|
============================= |
|
8
|
|
|
|
|
9
|
|
|
Walks through the steps of Natural Neighbor interpolation to validate that the algorithmic |
|
10
|
|
|
approach taken in MetPy is correct. |
|
11
|
|
|
""" |
|
12
|
|
|
########################################### |
|
13
|
|
|
# Find natural neighbors visual test |
|
14
|
|
|
# |
|
15
|
|
|
# A triangle is a natural neighbor for a point if the |
|
16
|
|
|
# `circumscribed circle <https://en.wikipedia.org/wiki/Circumscribed_circle>`_ of the |
|
17
|
|
|
# triangle contains that point. It is important that we correctly grab the correct triangles |
|
18
|
|
|
# for each point before proceeding with the interpolation. |
|
19
|
|
|
# |
|
20
|
|
|
# Algorithmically: |
|
21
|
|
|
# |
|
22
|
|
|
# 1. We place all of the grid points in a KDTree. These provide worst-case O(n) time |
|
23
|
|
|
# complexity for spatial searches. |
|
24
|
|
|
# |
|
25
|
|
|
# 2. We generate a `Delaunay Triangulation <https://docs.scipy.org/doc/scipy/ |
|
26
|
|
|
# reference/tutorial/spatial.html#delaunay-triangulations>`_ |
|
27
|
|
|
# using the locations of the provided observations. |
|
28
|
|
|
# |
|
29
|
|
|
# 3. For each triangle, we calculate its circumcenter and circumradius. Using |
|
30
|
|
|
# KDTree, we then assign each grid a triangle that has a circumcenter within a |
|
31
|
|
|
# circumradius of the grid's location. |
|
32
|
|
|
# |
|
33
|
|
|
# 4. The resulting dictionary uses the grid index as a key and a set of natural |
|
34
|
|
|
# neighbor triangles in the form of triangle codes from the Delaunay triangulation. |
|
35
|
|
|
# This dictionary is then iterated through to calculate interpolation values. |
|
36
|
|
|
# |
|
37
|
|
|
# 5. We then traverse the ordered natural neighbor edge vertices for a particular |
|
38
|
|
|
# grid cell in groups of 3 (n - 1, n, n + 1), and perform calculations to generate |
|
39
|
|
|
# proportional polygon areas. |
|
40
|
|
|
# |
|
41
|
|
|
# Circumcenter of (n - 1), n, grid_location |
|
42
|
|
|
# Circumcenter of (n + 1), n, grid_location |
|
43
|
|
|
# |
|
44
|
|
|
# Determine what existing circumcenters (ie, Delaunay circumcenters) are associated |
|
45
|
|
|
# with vertex n, and add those as polygon vertices. Calculate the area of this polygon. |
|
46
|
|
|
# |
|
47
|
|
|
# 6. Increment the current edges to be checked, i.e.: |
|
48
|
|
|
# n - 1 = n, n = n + 1, n + 1 = n + 2 |
|
49
|
|
|
# |
|
50
|
|
|
# 7. Repeat steps 5 & 6 until all of the edge combinations of 3 have been visited. |
|
51
|
|
|
# |
|
52
|
|
|
# 8. Repeat steps 4 through 7 for each grid cell. |
|
53
|
|
|
import matplotlib.pyplot as plt |
|
54
|
|
|
import numpy as np |
|
55
|
|
|
from scipy.spatial import ConvexHull, Delaunay, delaunay_plot_2d, Voronoi, voronoi_plot_2d |
|
56
|
|
|
from scipy.spatial.distance import euclidean |
|
57
|
|
|
|
|
58
|
|
|
from metpy.gridding import polygons, triangles |
|
59
|
|
|
from metpy.gridding.interpolation import nn_point |
|
60
|
|
|
|
|
61
|
|
|
plt.rcParams['figure.figsize'] = (15, 10) |
|
62
|
|
|
|
|
63
|
|
|
########################################### |
|
64
|
|
|
# For a test case, we generate 10 random points and observations, where the |
|
65
|
|
|
# observation values are just the x coordinate value times the y coordinate |
|
66
|
|
|
# value divided by 1000. |
|
67
|
|
|
# |
|
68
|
|
|
# We then create two test points (grid 0 & grid 1) at which we want to |
|
69
|
|
|
# estimate a value using natural neighbor interpolation. |
|
70
|
|
|
# |
|
71
|
|
|
# The locations of these observations are then used to generate a Delaunay triangulation. |
|
72
|
|
|
np.random.seed(100) |
|
73
|
|
|
|
|
74
|
|
|
pts = np.random.randint(0, 100, (10, 2)) |
|
75
|
|
|
xp = pts[:, 0] |
|
76
|
|
|
yp = pts[:, 1] |
|
77
|
|
|
z = (pts[:, 0] * pts[:, 0]) / 1000 |
|
78
|
|
|
|
|
79
|
|
|
tri = Delaunay(pts) |
|
80
|
|
|
delaunay_plot_2d(tri) |
|
81
|
|
|
|
|
82
|
|
|
for i, zval in enumerate(z): |
|
83
|
|
|
plt.annotate('{} F'.format(zval), xy=(pts[i, 0] + 2, pts[i, 1])) |
|
84
|
|
|
|
|
85
|
|
|
sim_gridx = [30., 60.] |
|
86
|
|
|
sim_gridy = [30., 60.] |
|
87
|
|
|
|
|
88
|
|
|
plt.plot(sim_gridx, sim_gridy, '+', markersize=10) |
|
89
|
|
|
plt.axes().set_aspect('equal', 'datalim') |
|
90
|
|
|
plt.title('Triangulation of observations and test grid cell ' |
|
91
|
|
|
'natural neighbor interpolation values') |
|
92
|
|
|
|
|
93
|
|
|
members, tri_info = triangles.find_natural_neighbors(tri, list(zip(sim_gridx, sim_gridy))) |
|
94
|
|
|
|
|
95
|
|
|
val = nn_point(xp, yp, z, (sim_gridx[0], sim_gridy[0]), tri, members[0], tri_info) |
|
96
|
|
|
plt.annotate('grid 0: {:.3f}'.format(val), xy=(sim_gridx[0] + 2, sim_gridy[0])) |
|
97
|
|
|
|
|
98
|
|
|
val = nn_point(xp, yp, z, (sim_gridx[1], sim_gridy[1]), tri, members[1], tri_info) |
|
99
|
|
|
plt.annotate('grid 1: {:.3f}'.format(val), xy=(sim_gridx[1] + 2, sim_gridy[1])) |
|
100
|
|
|
|
|
101
|
|
|
|
|
102
|
|
|
########################################### |
|
103
|
|
|
# Using the circumcenter and circumcircle radius information from |
|
104
|
|
|
# :func:`metpy.gridding.triangles.find_natural_neighbors`, we can visually |
|
105
|
|
|
# examine the results to see if they are correct. |
|
106
|
|
|
def draw_circle(x, y, r, m, label): |
|
107
|
|
|
nx = x + r * np.cos(np.deg2rad(list(range(360)))) |
|
108
|
|
|
ny = y + r * np.sin(np.deg2rad(list(range(360)))) |
|
109
|
|
|
plt.plot(nx, ny, m, label=label) |
|
110
|
|
|
|
|
111
|
|
|
|
|
112
|
|
|
members, tri_info = triangles.find_natural_neighbors(tri, list(zip(sim_gridx, sim_gridy))) |
|
113
|
|
|
delaunay_plot_2d(tri) |
|
114
|
|
|
plt.plot(sim_gridx, sim_gridy, 'ks', markersize=10) |
|
115
|
|
|
|
|
116
|
|
|
for i, info in tri_info.items(): |
|
117
|
|
|
x_t = info['cc'][0] |
|
118
|
|
|
y_t = info['cc'][1] |
|
119
|
|
|
|
|
120
|
|
|
if i in members[1] and i in members[0]: |
|
121
|
|
|
draw_circle(x_t, y_t, info['r'], 'm-', str(i) + ': grid 1 & 2') |
|
122
|
|
|
plt.annotate(str(i), xy=(x_t, y_t), fontsize=15) |
|
123
|
|
|
elif i in members[0]: |
|
124
|
|
|
draw_circle(x_t, y_t, info['r'], 'r-', str(i) + ': grid 0') |
|
125
|
|
|
plt.annotate(str(i), xy=(x_t, y_t), fontsize=15) |
|
126
|
|
|
elif i in members[1]: |
|
127
|
|
|
draw_circle(x_t, y_t, info['r'], 'b-', str(i) + ': grid 1') |
|
128
|
|
|
plt.annotate(str(i), xy=(x_t, y_t), fontsize=15) |
|
129
|
|
|
else: |
|
130
|
|
|
draw_circle(x_t, y_t, info['r'], 'k:', str(i) + ': no match') |
|
131
|
|
|
plt.annotate(str(i), xy=(x_t, y_t), fontsize=9) |
|
132
|
|
|
|
|
133
|
|
|
plt.axes().set_aspect('equal', 'datalim') |
|
134
|
|
|
plt.legend() |
|
135
|
|
|
|
|
136
|
|
|
########################################### |
|
137
|
|
|
# What?....the circle from triangle 8 looks pretty darn close. Why isn't |
|
138
|
|
|
# grid 0 included in that circle? |
|
139
|
|
|
x_t, y_t = tri_info[8]['cc'] |
|
140
|
|
|
r = tri_info[8]['r'] |
|
141
|
|
|
|
|
142
|
|
|
print('Distance between grid0 and Triangle 8 circumcenter:', |
|
143
|
|
|
euclidean([x_t, y_t], [sim_gridx[0], sim_gridy[0]])) |
|
144
|
|
|
print('Triangle 8 circumradius:', r) |
|
145
|
|
|
|
|
146
|
|
|
########################################### |
|
147
|
|
|
# Lets do a manual check of the above interpolation value for grid 0 (southernmost grid) |
|
148
|
|
|
# Grab the circumcenters and radii for natural neighbors |
|
149
|
|
|
cc = np.array([tri_info[m]['cc'] for m in members[0]]) |
|
150
|
|
|
r = np.array([tri_info[m]['r'] for m in members[0]]) |
|
151
|
|
|
|
|
152
|
|
|
print('circumcenters:\n', cc) |
|
153
|
|
|
print('radii\n', r) |
|
154
|
|
|
|
|
155
|
|
|
########################################### |
|
156
|
|
|
# Draw the natural neighbor triangles and their circumcenters. Also plot a `Voronoi diagram |
|
157
|
|
|
# <https://docs.scipy.org/doc/scipy/reference/tutorial/spatial.html#voronoi-diagrams>`_ |
|
158
|
|
|
# which serves as a complementary (but not necessary) |
|
159
|
|
|
# spatial data structure that we use here simply to show areal ratios. |
|
160
|
|
|
# Notice that the two natural neighbor triangle circumcenters are also vertices |
|
161
|
|
|
# in the Voronoi plot (green dots), and the observations are in the the polygons (blue dots). |
|
162
|
|
|
vor = Voronoi(list(zip(xp, yp))) |
|
163
|
|
|
voronoi_plot_2d(vor) |
|
164
|
|
|
|
|
165
|
|
|
nn_ind = np.array([0, 5, 7, 8]) |
|
166
|
|
|
z_0 = z[nn_ind] |
|
167
|
|
|
x_0 = xp[nn_ind] |
|
168
|
|
|
y_0 = yp[nn_ind] |
|
169
|
|
|
|
|
170
|
|
|
for i in range(len(nn_ind)): |
|
171
|
|
|
plt.annotate('{}, {}: {:.3f} F'.format(x_0[i], y_0[i], z_0[i]), xy=(x_0[i], y_0[i])) |
|
172
|
|
|
|
|
173
|
|
|
plt.plot(sim_gridx[0], sim_gridy[0], 'k+', markersize=10) |
|
174
|
|
|
plt.annotate('{}, {}'.format(sim_gridx[0], sim_gridy[0]), xy=(sim_gridx[0] + 2, sim_gridy[0])) |
|
175
|
|
|
plt.plot(cc[:, 0], cc[:, 1], 'ks', markersize=15, fillstyle='none', |
|
176
|
|
|
label='natural neighbor\ncircumcenters') |
|
177
|
|
|
|
|
178
|
|
|
for i in range(len(cc)): |
|
179
|
|
|
plt.annotate('{:.3f}, {:.3f}'.format(cc[i, 0], cc[i, 1]), xy=(cc[i, 0] + 1, cc[i, 1] + 1)) |
|
180
|
|
|
|
|
181
|
|
|
tris = tri.points[tri.simplices[members[0]]] |
|
182
|
|
|
for triangle in tris: |
|
183
|
|
|
x = [triangle[0, 0], triangle[1, 0], triangle[2, 0], triangle[0, 0]] |
|
184
|
|
|
y = [triangle[0, 1], triangle[1, 1], triangle[2, 1], triangle[0, 1]] |
|
185
|
|
|
plt.plot(x, y, ':', linewidth=2) |
|
186
|
|
|
|
|
187
|
|
|
plt.legend() |
|
188
|
|
|
plt.axes().set_aspect('equal', 'datalim') |
|
189
|
|
|
|
|
190
|
|
|
|
|
191
|
|
|
def draw_polygon_with_info(polygon, off_x=0, off_y=0): |
|
192
|
|
|
"""Draw one of the natural neighbor polygons with some information.""" |
|
193
|
|
|
pts = np.array([polygon[i] for i in ConvexHull(polygon).vertices]) |
|
194
|
|
|
for i in range(len(pts)): |
|
195
|
|
|
plt.plot([pts[i][0], pts[(i + 1) % len(pts)][0]], |
|
196
|
|
|
[pts[i][1], pts[(i + 1) % len(pts)][1]], 'k-') |
|
197
|
|
|
|
|
198
|
|
|
avex = np.mean(pts[:, 0]) |
|
199
|
|
|
avey = np.mean(pts[:, 1]) |
|
200
|
|
|
plt.annotate('area: {:.3f}'.format(polygons.area(pts)), xy=(avex + off_x, avey + off_y), |
|
201
|
|
|
fontsize=12) |
|
202
|
|
|
|
|
203
|
|
|
|
|
204
|
|
|
cc1 = triangles.circumcenter((53, 66), (15, 60), (30, 30)) |
|
205
|
|
|
cc2 = triangles.circumcenter((34, 24), (53, 66), (30, 30)) |
|
206
|
|
|
draw_polygon_with_info([cc[0], cc1, cc2]) |
|
207
|
|
|
|
|
208
|
|
|
cc1 = triangles.circumcenter((53, 66), (15, 60), (30, 30)) |
|
209
|
|
|
cc2 = triangles.circumcenter((15, 60), (8, 24), (30, 30)) |
|
210
|
|
|
draw_polygon_with_info([cc[0], cc[1], cc1, cc2], off_x=-9, off_y=3) |
|
211
|
|
|
|
|
212
|
|
|
cc1 = triangles.circumcenter((8, 24), (34, 24), (30, 30)) |
|
213
|
|
|
cc2 = triangles.circumcenter((15, 60), (8, 24), (30, 30)) |
|
214
|
|
|
draw_polygon_with_info([cc[1], cc1, cc2], off_x=-15) |
|
215
|
|
|
|
|
216
|
|
|
cc1 = triangles.circumcenter((8, 24), (34, 24), (30, 30)) |
|
217
|
|
|
cc2 = triangles.circumcenter((34, 24), (53, 66), (30, 30)) |
|
218
|
|
|
draw_polygon_with_info([cc[0], cc[1], cc1, cc2]) |
|
219
|
|
|
|
|
220
|
|
|
########################################### |
|
221
|
|
|
# Put all of the generated polygon areas and their affiliated values in arrays. |
|
222
|
|
|
# Calculate the total area of all of the generated polygons. |
|
223
|
|
|
areas = np.array([60.434, 448.296, 25.916, 70.647]) |
|
224
|
|
|
values = np.array([0.064, 1.156, 2.809, 0.225]) |
|
225
|
|
|
total_area = np.sum(areas) |
|
226
|
|
|
print(total_area) |
|
227
|
|
|
|
|
228
|
|
|
########################################### |
|
229
|
|
|
# For each polygon area, calculate its percent of total area. |
|
230
|
|
|
proportions = areas / total_area |
|
231
|
|
|
print(proportions) |
|
232
|
|
|
|
|
233
|
|
|
########################################### |
|
234
|
|
|
# Multiply the percent of total area by the respective values. |
|
235
|
|
|
contributions = proportions * values |
|
236
|
|
|
print(contributions) |
|
237
|
|
|
|
|
238
|
|
|
########################################### |
|
239
|
|
|
# The sum of this array is the interpolation value! |
|
240
|
|
|
interpolation_value = np.sum(contributions) |
|
241
|
|
|
function_output = nn_point(xp, yp, z, (sim_gridx[0], sim_gridy[0]), tri, members[0], tri_info) |
|
242
|
|
|
|
|
243
|
|
|
print(interpolation_value, function_output) |
|
244
|
|
|
|
|
245
|
|
|
########################################### |
|
246
|
|
|
# The values are slightly different due to truncating the area values in |
|
247
|
|
|
# the above visual example to the 3rd decimal place. |
|
248
|
|
|
|
Unidata /
MetPy
